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Preface 


This  study  shows  that  optimization  techniques  can  be 
applied  to  the  solution  of  complex,  boundary-layer  flow  prob- 
lems.  It  builds  on  previous  work  done  by  students  and  facul¬ 
ty  at  the  Air  Force  Institute  of  Technology.  Captain  Karen 
Lange  working  with  Major  James  K.  Hodge  developed  a  computer 
code  to  solve  boundary- layer  problema.  Captain  Lange'a  com¬ 
puter  code  was  modified,  and  optimization  codes  studied  by 
1st  Lieutenant  Bruce  K.  Boyd  and  Captain  Sal  A.  Leone  were 
added.  Without  the  work  of  these  individuals,  my  extension 
of  both  studies  to  optimization  of  two-dimensional,  boundary- 
layer  code  would  not  have  been  possible. 

I  would  like  to  take  this  opportunity  to  especially 
thank  my  aidvisor.  Dr.  Sal  A.  Leone,  for  his  guidance,  know¬ 
ledge,  and  assistance  throughout  this  study.  Of  course, 
this  study  could  not  have  proceeded  very  far,  without  the 
knowledge  and  expertise  of  Major  Hodge  in  Computer  Fluid 
Dynamics.  I  thank  him  for  always  willingly  sharing  his  ex¬ 
pertise  when  I  needed  it.  Lastly,  I  thank  Lt.  Col.  Eric 
Jumper  and  Dr.  M.  L.  Rasmussen  for  giving  me  the  foundation 
in  theoretical,  fluid  dynamics  that  enabled  me  to  understand 
the  physics  of  the  flow  problem.  If  in  some  small  way,  this 
study  helps  further  studies  of  similar  problems  or  advances 
the  study  of  hypersonic,  fluid  dynamics,  it  is  a  success. 

Donald  E.  Coffey,  Jr. 
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This  study  develops  sn  sdapt  ive-grid  method  which  mini 
mixes  the  truncation  error  in  the  finite-difference  solution 
The  study  solves  compress  ibie ,  steady-state,  boundary- layer 
equations  assuming  perfect-gas  flow  over  an  isothermal  wall. 
The  Dorodnitsyn,  compressibility  transformation  changes  the 
bounds ry- 1  aye r  equations,  as  expressed  in  two-dimensional, 
cartesian  coordinates,  into  an  incompressible  form.  The 
equations  are  then  transformed  into  variables  of  a  compu¬ 
tational  plane.  Implicit  Sue c e s * ive-Over-Be 1  ax a t ion  (SOB) 
solves  the  f in i t e-d i f f e r enc ed .  computational,  bounds ry- 1  ay e r 
equations.  Comparison  of  the  computed  solution  for  incom¬ 
pressible  flow  over  a  flat  plate  to  Blasius',  exact  solution 
shows  the  boundary- layer  code  is  accurate. 

The  adaptive-grid  method  uses  Powell's  method  to  opti¬ 
mize  the  solution  grid  by  minimizing  the  sum  of  the  third 
derivative  in  the  computational  plane  of  the  tangential  velo¬ 
city  component.  Powell's  method  finds  the  grid,  control 

function,  Q,  in  an  elliptic,  grid  equation,  Y _  +  QY  =0 , 

nn  n 

which  minimizes  a  specified  function.  The  grid  equation 
generates  the  grid  spacing  at  the  end  of  the  plate.  This 
spacing  is  then  streamwise  scaled  across  the  remaining  grid. 
Minimizing  the  sum  of  the  square  of  the  third  derivative  in 
the  computational  plane  of  the  tangential  velocity  component 
,  over  the  entire  domain  decreases  the  truncation  error 


the  beet  of  the  functions  tested.  This  stndy  tests  the  sun 
of  the  squares  of  the  first,  second,  and  third  derivative  of 
the  tangential  velocity  as  minimized  functions.  The  accuracy 
of  the  optimized,  adaptive-grid  solution  is  greater  than  the 
original,  fixed-grid  solution.  The  study  then  applies  the 
optimization  to  supersonic  and  hypersonic  problems.  The 
computed,  adaptive-grid  solutions  show  good  correlation  with 
theoretical  models  for  supersonic  and  hypersonic  flow  devel¬ 
oped  by  Van  Driest  and  Eckert. 


ADAPTIVE-GRID  OPTIMIZATION  FOR 


MINIMIZING  STEADY-STATE,  TRUNCATION  ERROR 


Chapter  Ij 


Every  flight  of  the  Space  Shuttle  highlights  the  know¬ 
ledge  aeronautical  engineers  have  of  hypersonic  flight.  The 
success  of  each  flight  shows  a  basic  understanding  of  the 
problems  involved  with  flight  at  high  altitudes  and  at  speeds 
above  Mach  S.  However,  if  future  aeronautical  engineers  are 
to  build  better  shuttles  or  to  develop  a  transatmospher ic 
vehicle,  increased  understanding  of  hypersonic  flight  will  be 
necessary.  Presently,  experimentation  in  the  wind  tunnel  and 
Space  Shuttle  flights  validates  portions  of  the  hypersonic 
theory.  Unfortunately,  experimental  work  of  this  sort  is 
very  expensive.  Newer  methods  of  modeling  flow  using  compu¬ 
ters  offer  hope  for  simulating  characteristic,  fluid  proper¬ 
ties  in  the  hypersonic  regime  without  the  experimental  ex¬ 
pense.  Using  computational  fluid  methods,  unexpected  exper¬ 
imental  results  can  be  confirmed,  and,  conversely,  experi¬ 
ments  can  be  designed  to  validate  results  from  computational 
fluid  modeling.  In  this  way,  better  understanding  of  hyper¬ 
sonic  fluid  behavior  can  be  gained  and  can  lead  to  more  effi¬ 
cient  hypersonic  vehicles. 

An  example  of  the  usefulness  of  computational  modeling 
for  explaining  the  results  of  experimental  data  is  the  non- 
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isothermal  vail  effects  which  have  occurred  on  the  Space 
Shnttle.  The  surface  of  any  body  travelling  at  high  speeds 
through  any  fluid  heats  up.  The  effect  the  heated  surface 
has  on  the  flov  field  or  boundary  layer  around  the  vehicle 
depends  on  the  heat  transfer  coefficient  of  the  surface.  If 
the  body  is  one  material  and  is  heated  evenly,  the  convection 
which  takes  place  is , isothermal .  However,  if,  as  is  the  case 
with  the  Space  Shuttle,  there  are  several  different  materials 
joined  together  on  the  same  surface,  non- isothermal  convec¬ 
tion  o  jura.  Non- isothermal  convection  can  change  the  air 
and  heat  flow  on  the  vehicle.  Boundary- layer  theory  for  the 
hypersonic,  non- iso  the rma 1  wall  predicts  a  discontinuity  in 
wall  temperature  in  the  area  of  a  material  change  (11).  How¬ 
ever,  wind  tnnnel  data  suggests  a  much  slower,  wall-temper¬ 
ature  recovery  than  theory  suggests  (21:2).  Roberts  and 
Lange,  in  two  separate  studies,  investigated  the  non-isother- 
■al-wall  effect  using  computational  fluid  dynamics  (CFD) 
(15,21).  Roberts  modeled  flows  with  two-dimensional  Navier- 
Stokes  equations.  Lange  modeled  flow  using  an  unsteady, 
boundary- 1 ayer  analysis.  Both  studies  showed  a  non-isother- 
mal,  wall  effect.  However,  in  the  case  of  the  boundary-layer 
analysis,  inaccuracies  in  the  method  tended  to  smear  the 
quantitative  results.  The  disadvantage  with  most  computer, 
modeling  techniques  is  truncation  error  and  round-off  error 
smear  the  results  in  h i gh- g r ad ien t  areas.  These  high-grad- 
ient  areas  also  are  the  parts  of  the  flov  where  the  non-iso- 
thermal,  wall  effect  takes  place..  Therefore,  the  engineer 
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■nit  decrease  the  errors  is  the  high-gradient  areas  to  aore 
accurately  model  flow  conditions.  In  this  way,  computer 
fluid  modeling  can  more  accurately  predict  flow  conditions 
such  as  the  non- isothermal ,  wall  effect,  which  otherwise 
would  not  be  seen. 

Background 

Using  an  adaptive  grid  which  is  changed  after  each  solu¬ 
tion  of  the  modeling  equations  is  a  way  to  increase  the  ac¬ 
curacy  of  the  computer  solution  without  increasing  the  number 
of  solution,  grid  points.  The  easiest  way  to  increase  reso¬ 
lution  of  regions  where  fluid  properties,  like  temperature 
and  pressure,  are  changing  rapidly  is  to  increase  the  number 
of  grid  points.  The  other  way  is  to  use  an  adaptive  grid  for 
solving  the  model  equations.  The  adaptive  grid  concentrates 
grid  points  in  h igh-g r ad i en t  areas,  (i.e.  boundary  layers), 
to  increase  resolution  and  spreads  out  grid  points  in  low- 
gradient  regions.  This  keeps  the  total  number  of  grid  points 
small,  but  placos  grid  points  to  produce  the  best  resolu¬ 
tion.  For  a  time-dependent  problem,  Ghia,  Ghia,  and  Shin 
develop  a  flow-dependent,  adap t ive- gr id  method  to  increase 
the  accuracy  of  the  computer  solution.  This  adaptive-grid 
method  reforms  the  solution  grid  after  each  time  iteration  of 
the  modeling  equations  by  minimizing  the  magnitude  of  the 
convective  terms  of  the  governing  equations  which  indirectly 
minimizes  the  truncation  error  in  the  governing  equations 
(7:36-37).  Their  treatment  of  the  two-dimensional,  boundary- 
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layer  equations  adapts  the  grid  in  normal  and  streanwise 
directions.  However,  since  the  gradients  in  a  boundary  layer 
are  not  very  large  in  the  streanwise  direction,  it  may  be 
possible  to  optimize  the  solution  grid  in  only  the  n.raal 
direction  and  still  get  an  accurate  solution.  If  the  parti¬ 
cular,  flow  solutions  do  not  have  similar,  velocity  profiles, 
the  solution  grid  can  be  coupled  with  a  parabolic-grid  genet- 
ition  method  (9).  Regardless,  one-dimensional,  grid,  optimi¬ 
zation  techniques  can  be  applied  to  optimize  the  grid  in  the 
normal  direction.  Leone  and  Hodge  suggest  optimizing  the 
adaptive  grid  in  one  dimension  with  Powell's  minimization 
method  (16).  Boyd  also  uses  a  method  to  optimize  the  adap¬ 
tive  grid  in  one  dimension.  Using  a  least  squares  curve  fit 
to  model  a  quadratic  equation  and  a  Newton-Raphson  optimiza¬ 
tion  of  a  grid  control  function,  Q,  Boyd  attempts  to  minimize 
the  tru-cation  error  of  the  third  derivative  of  the  dependent 
vttisbi  ■  to  get  a  more  accurate  one-dimensional,  flow  solu¬ 
tion  (3:11).  After  optimizing  the  solution  grid  in  the 
normal  direction  with  one-dimensional  techniques,  the  charac¬ 
teristics  of  the  boundary  layer  solution  can  be  used  to  scale 
the  normal  optimization  in  the  streamwise  direction.  The 
resulting  optimized,  solution  grid  generates  more  accurate 
flow  solutions. 

Objective 

This  thesis  develops  a  method  for  optimizing,  a  steady- 
state,  adap t ive-gr id  solution  of  two-dimensional,  flow  prob- 
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lens  by  minimizing  the  truncation  error  in  the  streamwise 
velocity  component.  The  adaptive  grid  minimizes  the  snm  of 
t.he  square  of  the  third  derivative  of  the  streamwise  velocity 
in  the  transformed  plane,  n*in*  Powell's  minimization 

method.  Powell's  method  iterates  on  an  n-d imens ional  array 
using  conjugate  directions  until  finding  the  minimum  of  the 
function  without  calculating  the  derivatives  of  the  function 
on  that  array.  In  this  study,  Powell's  method  primarily 
minimizes  the  third  derivative  of  the  streamwise  velocity 
component  by  optimizing  the  grid  control  function,  Q,  in  the 
equation,  Y^  +  QY^  =0 .  Because  of  the  flexibility  of  Pow¬ 
ell's  method,  the  effect  of  optimizing  the  solution  grid  by 
minimizing  the  first,  second,  or  a  combination  of  all  three 
streamwise  velocity  derivatives  is  checked  to  verify  trunca¬ 
tion  error  arguments.  Regardless  of  the  specific  minimized 
function  used,  the  method  optimizes  the  grid  in  the  normal 
direction.  The  grid  is  then  scaled  in  the  streamwise  direc¬ 
tion  using  boundary- 1  ay e r  theory.  A  modification  of  Lange's 
boundary- layer  code  then  solves  the  flow  problem.  Comparing 
the  incompressible  and  compressible  flow  solutions  to  exact 
solutions  and  Lange's  non-adaptive  numerical  solutions  veri¬ 
fies  the  accuracy  of  the  adaptive-grid  solution.  The  in¬ 
crease  in  accuracy  with  the  optimized,  adaptive  grid  should 
be  evident  in  the  velocity  calculations  and  in  the  comparison 
of  heat  transfer  along  the  surface. 
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Chanter  I I :  Theory 


Boundary-Laver  Eouat ions 

The  complete  Nayier-Stokee  equations  predict  fluid  char¬ 
acteristics  for  all  flow  conditions  (23:64).  This  set  of 
equations  is  however  very  difficult  to  solve;  Since  the 
majority  of  viscous  effects  and  heat  transfer  effects  appear 
in  a  thin  layer  adjacent  to  the  flow  surface,  this  study 
investigates  only  the  portion  of  flow  in  that  region,  known 
as  the  boundary  layer.  Flow  properties  in  the  boundary  layer 
allow  several  simplifying  assumptions  which  result  in  a  sim¬ 
pler  set  of  flow,  modeling  Aquations  than  the  Navier-Stokes 
equations.  For  high  Reynold's  number  flows,  viscous  effects 
in  the  form  of  shearing  itresi  at  the  wall,  rw  »  p ( d u/ dy ) y_0 , 
and  velocity  gradients  in  the  normal  direction,  du/dy,  are 
very  large  inside  the  boundary  layer.  Outside  the  boundary 
layer,  flow  is  essentially  inviscid  with  negligible,  velocity 
gradients  (23:128-129).  .  Therefore,  within  the  boundary  layer 
thickness,  6,  it  is  assumed  that: 

a.  gradients  in  the  normsl  directions  are  much  larger 
than  gradients  in  the  streamwise  direction, 

b. the  normal  velocity,  v,  is  much  smaller  than  the 
streamwise  velocity,  u,  and 

c.  all  terms  of  order  S/L  or  smaller  are  negligible. 

For  laminar  flow  over  a  flat  plate  or  wedge,  the  pressure 
gradient  in  the  streamwise  direction,  dp/dx,,  of  the  iuviscid 


region  is  impressed  on  the  boundary  layer.  This  pressure 
gradient  is  assumed  to  be  zero  for  a  flat  plate  parallel  to 
the  flow.  Further,  the  boundary  layer  equations  are  non-di¬ 
mens  ional  ized  using  the  following  relationships. 


x“  xi  ,  y-  zl>  t*  _ t ' D 1 _  .  u-  aJ _ ,  ▼  -  v  *  , 

L'  L'  L'  U*  U' 


H=  H  *  T=>  Co  *  T  *  p«  o  1 

UV  OV  p«a  '  D  1 


»»"  ilI _ »  pm  nl 

Pm'  Pm 


The  resulting  non-dimens ional ized ,  b ound ary- 1  ay e r  equations 
for  compressible,  laminar  flow  are: 


Continuity: 


3  o  +  3  ou  +  3 ov 
3 1  3x  3 y 


Momentum: 


3ou  +  3  ouz  +  3  ouv 
3 1  3x  3y 


-is.  +  il_  is 

3x  3y  Re  3y 


Energy : 


3  o  »  0 
3y 

3pH+3puH+3pvH  *  3  p  +  3  u  3  H 
3t  3x  3y  3t  3y  RePr  3y 


♦  u(Pr-l)  3u  12  (4) 

3y  RePr  3y 

Boundary  conditions  for  the  boundary  layer  equations  are 
u'=  v ' =  0  at  y ’=  0  and  u ' =  D#  at  y'=6.  Assuming  incompressi¬ 
ble  flow  or  transforming  the  compressible  equations  into  an 
incompressible  form  using  the  Dorodnitsyn  transformation 
further  simplifies  the  above  equations  (22:101). 


ip r e  s  s  ib 1 e .  Boundary-Layer  Theory 

If  the  flow  density  is  constant,  the  velocity  profile. 
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as  well  as  the  heat  transfer  characteristics  near  the  wall, 
have  already  been  derived  theoretically  for  both  exact  and 
approximate  cases.  In  the  incompressible  case,  the  boundary 
layer  equations  are: 


Cont inuity : 

O+O  -  o 

31  3  Y 

(5) 

Momentum: 

4JJ+UO+V32  = 

-1  3i+3  I 

3U\ 

(6a) 

3 1  3X  3  Y 

p  3X  3Y  ' 

^  Re 

3Y  / 

it  «  0 

3  Y 

■  • 

(6b) 

Energy:  3H+U3JI+V3H  3  L  3n  +  3  f  uf)  3  H  | 

3 1  31  3 Y  p  3 1  3Y\RePr  3Y / 

+  3  /  pu.  ( Pr-1 )  3P*/2  \  (7) 

3Y\  RePr  3Y  J 

where  u=U,  v-V,  x  *X,  and  y“ Y 

Blasins  developed  the  exact  solution  to  incompressible,  boun¬ 
dary  layer  flow  over  a  flat  plate.  The  results  of  hia  exact 
solution  define  the  thickness  of  the  boundary  layer.  If  the 
edge  of  the  boundary  layer  is  assumed  to  be  where  U’/U.  is 
.994,  the  boundary  layer  thickness,  3,  at  any  streamwise 
location,  x',  is 


5(x')  “  5 ■ 2  x  f .  (8) 

where  Rex  =  pa>  *  u  *  x  *  /  p  ’  m .  Blasins'  solution  also  gives  an 
exact  value  for  the  skin  friction  coefficient,  (?£. 


(9) 
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The  computer  solution  of  the  incompressible,  fist,  piste 
problem  should  mstch  this  exact  Cj  result. 

Though  Blasius*  incompressible,  flat,  plate  solution  is 
exact  and  is  the  best  description  of  this  velocity  profile, 
the  von  Zarman-Pohlhausen,  approximate  solution  of  the  two- 
dimensional,  flow  problem  offers  more  easily  obtainable  velo¬ 
city  profile  results.  The  von  Karman-Pohlhausen  method 
solves  the  momentum- integral  equation  (23:160).  The  velocity 
distribution  is  a  fourth-order  polynomial  where  w-y  *  /  6  ( x)  . 

0  -  01  *  aw  +  b*»^  +  +  din*  (10) 

The  four  boundary  conditions  required  to  solve  the  coeffi¬ 
cients  in  Eq  (10)  are 

y'-  0  :  0  -  o  a£o  -  1  a*  -  -o#  40, 

ar7  p  ax  dx 

y  *  *  8:0-1  dO  -  0,  3*0  -  0  (11) 

ai  aT* 

The  coefficients,  a-d,  in  Eq  (10)  are 

a-  2+0/6  ,  b*=  -0/2  ,  c=  -2+0/2,  d=  1-0/6 

where  0-  6*udO0  , (12) 

pdX 

(23:207).  0  is  the  Pohlhausen  pressure,  gradient  parameter. 

For  flat  plate  and  wedge  flows  with  no  streamwise,  pressure 
gradient,  0  is  zero.  The  von  Karman-Pohlhausen  solution 
provides  another  verification  for  numerically-derived,  velo¬ 
city  profiles  in  this  study.  After  verifying  the  accuracy  of 
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the  boundary- laye r  and  adaptive-grid  code  with  the  Blasius  or 
von  Karman-Pohlhausen  resalts,  then  compressible  forms  of  the 
boundary-layer  equations  are  solved. 

Compressibility  Considerations 

Since  the  incompressible,  boundary^ 1  ay e r  equations,  Eqs 
(5),  (6),  and  (7),  already  have  approximate  and  exact  solu¬ 
tions,  the  solution  of  the  compressible,  boundary- layer  equa¬ 
tions,  Eqs  (2),  (3),  and  (4),  is  possible  if  the  change  in 
density  and  viscosity  across  the  boundary  layer  can  be  ac¬ 
counted  for.  The  Dorodnitsyn  transformation  changes  the 
compressible  equations  into  an  incompressible  form.  The 
compressible  equations  can  then  be  solved  sis  incompressible 
equations.  The  Dorodnitsyn  transformation  is  *a  nonlinear 
stretching  of  the  normal  coordinate"  (20).  The  transformation 
takes  the  normal  coordinate,  y,  in  the  physical  space  and 
stretches  it  into  a  new  normal  coordinate,  T,  in  the  trans¬ 
formed  space.  The  transformed,  space  coordinates  are  now 
expressed  in  X  and  7  where 

X  -  x,  I  »  |  _jlL  dy  (13) 

JO  pe' 

The  transformed,  grid  coordinates  lead  to  a  change  in  the 
velocity  components,  as  shown  in  Append!  A.  The  transform¬ 
ed.  velocity  components  are 

U  -  u,  V  =  pv  +  3T/at  +  u3Y/3x  (14) 


Using  the  transformed  variables  in  Eqs  (13)  and  (14),  changes 


the  compressible,  boundary- 1  aye r  equations,  Eqs  (2)  -  (4), 
into  the  incompressible  form  in  Eqs  (5)  -  (7).  The  boundary 
layer  thickness,  S(x),  also  changes  into  a  transformed,  boun¬ 
dary-layer  thickness,  A(X).  Assuming  a  calorically-perf ect 
gas,  the  density  ratio,  p  '  /  p  '  a ,  equals  the  inverse  of  the 
specific  enthalpy  ratio,  hs  ’/hs*,  and 
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(15) 


An  explicit  expression  for  A(x)  comes  from  an  integral  analy¬ 
sis  of  the  new,  boundary-layer  equations  (20). 


where  B  =  dD/P 

d(Y/A) 

and  fQ  -  P/De(l-0/De)  d(Y/A) 

a,  in  the  above  equation,  is  0  for  two-dimensional  flow  and  1 
for  axisymmetric  flow  (Appendix  P) .  Assuming  flow  over  the 
flat  plate  or  wedge  develops  a  cubic,  velocity  profile  such 
that 

0/Be  =  1.5(Y/5)  -  0.5<Y/6)3  ,  (17) 

the  value  of  2  B / f q  is  280/13  .  Known  flow  conditions  then 
determine  A(X).  After  solving  the  flow  for  the  incompres¬ 
sible  form,  inveise  transformations  compute  the  solution  in 


the  physical  space.  The  inverse  t r ans f o rma t ions  result  from 
integrating  the  computed,  enthalpy  ratio  to  recover  the  phy¬ 
sical-grid,  normal  coordinate,  or 


y  “f  hs'  dT  (18) 

Jo  k*  e ' 


The  compressible  solution  of  the  bonrdary  layer  equations  in 
the  physical  plane  for  a  given  set  of  flow  conditions  is  then 
complete . 


Flow  Prcoerties 

Ob  I i a  n  e .  Shock-Wave  Theory. 

For  supersonic  flow  over  a  flat  plate,  the  flow  condi¬ 
tions  snch  as  density,  pressure,  temperature,  and  velocity 
remain  relatively  constant  ontside  the  boundary  layer.  How¬ 
ever.  If  the  flat  plate  is  inclined  to  the  direction  of  flow 
to  simulate  a  wedge,  the  flow  conditions  change  as  flow  pas¬ 
ses  through  the  shock  wave  which  develops  in  front  of  the 
wedge.  This  study  assumes  a  perfect}  gas  flow  with  no  dis¬ 
sociation  of  the  gas  molecules.  For  the  perfect  gas, 

p*  «  p  '  R  T' ,  and  p  *  v  2.  (19) 

y-1  T 

where  the  gas  constant,  K,  is  1715  f t- lb / s 1 ug-°R  and  the 
specific  heat  ratio,  y,  is  1.4.  Oblique,  shock-wave  theory 
applies  to  this  flow  problem.  To  find  the  conditions  behind 
the  shock  wave,  the  shock  angle,  0,  must  be  found.  Liepmann 
and  Roshko  use  an  explicit  form  of  the  shock-angle  equation 
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> V7 V  v.'i-v.v .  v.  v. ».* 


(17:87).  However,  the  implicit  form. 


M#2  sin2p  -  1 


-X±J^M  = 
2 


s inB  s in6 
co  * ( p-Q) 


(20) 


is  more  useful  for  iteratively  competing  the  shock  angle,  p, 
given  a  deflection  angle,  0,  and  a  freestream,  Mach  number, 
M.,.  With  the  shock  angle  and  freestream  values  of  Mach,  Mm, 
pressure,  pa,  and  temperature,  Tm,  the  equations 


M  2  s in2 ( p-6)  »  2  ±  (y-1)  sinzB 

‘  2M00*sinzp  -  (y-1) 

P.  *  P«[l  +  -lx  (M00*sin2p-1)  ) 


T  =  T 
e  < 


T  +  1 

1  4  «  2  R  -  < 


,fl  ♦  2  (v-1)  M_*  sin2B  .-(.tM_‘  s  in2P  ±1)1 
[  (r  +  1)  M.1  » in2  P  J 


(21) 

(22) 

(23) 


define  the  edge  conditions  behind  the  shock  wave  for  Mach 
number,  pressure,  and  temperature  (17:86).  All  other  edge 
conditions  such  as  enthalpy  and  density  can  be  found  from 
these  quantities.  These  edge  conditions  affect  the  solution 
of  the  boundary- layer  equations  over  the  given  surface. 


Viscosity. 

Viscosity  is  not  a  constant  across  the  boundary  layer. 
Since  viscosity  is  a  function  of  the  fluid  temperature  which 
changes  across  the  boundary  layer,  viscosity  also  changes. 
Generally,  engineers  use  Sutherland's,  viscosity  law. 


|i'  “  2 .2685  x  10~8  T'3/2  slues 
T'  +  198.6  °R  ft-sec 


(24) 


for  calculating  viscosity  as  the  fluid  temperature  changes, 
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However,  Fiore  suggests  Sutherland's  law  is  not  valid  for  the 
low  temperatures  resulting  from  hypersonic,  flow  experiments 
(6:56).  Flow  temperatures  in  hypersonic,  wind-tunnel  testing 
generally  range  from  30°R  to  200°R  (6:56).  A  more,  exact 
expression  for  low-temperature  viscosity  applies.  Keyes 
obtains  bistter  low-temperature  calculations  of  viscosity 
using  the  relation 

H  '  «  2.32  x  IQ"8  T'1'2  sluts  (25) 

1+I220/T'  ( 10* /  1  )  ]  ft-aec 

(13).  In  addition,  this  1 ow- t empe r a tur e ,  viscosity  law  gives 
virtually  the  same  results  for  temperatures  cbove  300°R. 
Cappelano  states  "variation  between  the  alternate  form  and 
Sutherland's  law  is  negligible  above  a  temperature  of  T» 

300°R  "  (4:20).  Below  300°R,  there  is  a  difference.  The 
new,  viscosity  law  is  more  accurate  for  low  temperatures. 
Therefore,  unless  computational  results  must  be  compared  with 
previous  studies  which  use  Sutherland's  law,  this  study  uses 
Keyes'  viscosity  law  where  appropriate. 

Convective  Heat  Transfer  Coefficient. 

Convection  is  the  primary  method  of  heat  transfer  across 
the  boundary  layer.  For  a  constant  temperature  surface,  the 
heat  flow,  q,  is  calculated  by  the  convective,  heat  relation 

q  «  hA(Tw  -  T#w)  (26) 

where  Tw  is  the  wall  temperature,  T#w  is  the  adiabatic,  wall 


temperature,  and  A  is  the  wall  surface  area.  Thus,  it  is  not 


£*«  A.  .A.  *  A 


K; 

k  s 

?  1 

i  * 


.  /*  /•  .**  ;*»  .A  vAi  .A  t*V  ^  A  ;•%  ^  * 


surprising  that  the  ability  to  accurately  compute  the  convec¬ 
tive,  heat,  transfer  coefficient,  h,  is  important  to  any 
boundary- layer  study.  Numerically,  the  computer  code  calcu¬ 
lates  the  heat  flow  rate  per  unit  area  from  the  flow  solution 


us  ing 


q/A  -  -k  3T’ /9y * 


where  k  is  the  thermal,  conductivity  coefficient  for  a  fluid. 
The  wall  temperature,  Tv  is  specified.  Taw  varies  with  Mach 
number  and  fluid  temperature  by  the  adiabatic,  wall  tempera¬ 


ture  , 


Te[l+(Pr)1/2  (v-1)M_21 
2 


The  Prandtl  number  is  assumed  constant  and  is 


Pr  *  Cp|i/k 


The  factor,  Pr  -  ,  in  Eq  (28)  is  the  adiabatic,  recovery 

factor  for  laminar  flows  (12:213).  In  assuming  an  adiabatic 
wall,  Schetz  explains, 

the  temperature  that  the  wall  attains  at  equili¬ 
brium  will  depend  on  how  much  of  this  kinetic  ener¬ 
gy  (kinetic  energy  of  the  flow]  is  recovered  on 
the  wall.  This  is  expressed  in  the  recovery  fac¬ 
tor...  (22:5) 

Solving  Eq  (26),  using  Eqs  (27)-(29),  gives  a  numerical  solu¬ 
tion  for  h.  There  is  also  a  theoretical  solution  for  h. 

The  theoretical  solution  comes  from  Blasius'  solution  of  the 
flat-plate  problem  as  well  as  Eckert's,  flat-plate  theory. 
From  Blasius',  exact  solution  for  flat  plate  flow,  Eq  (9) 
gives  a  value  for  Cf.  Cf  is  also  related  to  the  Stanton 
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number,  S  t  x ,  by 


fx  -  2StxPr2/3 

(30) 

St  -  h  * 

(31) 

‘  ~<fp  U. 

(12:195).  Substituting  Eq  (31)  into  Eq  (30),  equating  the 
result  to  Eq  (9),  and  solving  for  h  gives  a  theoretical  value 
of  h: 

V  *  *332  [|ie'Pe'V/*']1/2  (32) 

Pr* 

This  value  of  h  is  only  for  flows  where  viscosity  and  density 
are  constant,  the  wall  is  isothermal,  and  the  flnid  is  a 
perfect,  s ing 1 e- sp ec i e s  gas  with  no  dissociation  of  molec¬ 
ules.  It  applies  therefore  to  incompressible  flows  treated 
in  this  study.  However,  for  compressible,  boundary-layer 
flows,  where  density  and  viscosity  vary,  the  theory  must  be 
modified  slightly.  Eckert  suggests  the  use  of  a  reference 
temperature  which  is  representative  of  the  temperatures  ac¬ 
ross  the  boundary  layer  in  the  previous  constant  property 
relations  (22:96).  Defined  as 

T* '  -  Te'+.5(Tw'-T#')+.22(Taw>-Te')  (33) 

the  reference  temperature  is  used  to  calculate  reference 
values  of  viscosity  and  density,  p*'  and  p*',  respectively. 
Then,  Stanton  number  and  h  for  compressible  flows  are  calcu¬ 
lated  using  the  reference,  *,  conditions.  Another  reference 


16 


ia 


enthalpy,  hrejf,  it  calculated  by  using  pm,  pa,  end  ua  in 
place  of  the  edge  conditions  in  Eq  (32).  To  evaluate  how  h 
changes  fro*  its  freeitreaa  value,  the  theoretical  h  values 
for  the  incompressible  and  compressible  cases  are  compared  to 
the  freestreao  h.  For  incompressible  flow  over  a  flat 
plate,  with  edge  conditions  equal  to  freestream  conditions, 
h/hre£  should  be  one  for  both  numerical  and  theoretical  val¬ 
ues  of  h.  Bow  close  the  numerical  h/hre£  is  to  the  theoreti¬ 
cal  b/href  indicates  how  well  the  computed  solution  models 
the  exact,  incompressible  solution. 

Theoret ical  Limitat ions 

The  assumptions  which  form  the  theoretical  basis  for  the 
boundary-layer  theory  of  this  study  pose  several  limitations 
on  the  general  application  of  the  solutions  for  all  flow 
conditions  in  an  experimental  setting.  The  flows  are  as¬ 
sumed  to  be  perfect  gases.  This  limits  the  results  to  air 
flows  whose  temperatures  and  pressures  guarantee  negligible, 
rarefied-gas  effects,  i.e.  molecular  vibration  or  dissocia¬ 
tion.  Also,  the  requirement  of  a  zero  streamwise  pressure 
gradient  across  the  boundary  layer  limits  the  shapes  to  which 
it  can  be  applied  without  modification.  These  shapes  include 
flat  plates,  wedges,  and  cones.  Most  important  is  a  limita¬ 
tion  on  the  physical  effects  which  can  be  predicted  in  the 
hypersonic  regime.  Boundary- layer  theory  is  incapable  of 
predicting  the  effects  of  shock-boundary-layer  interaction. 
The  shock  wave  produced  at  the  leading  edge  of  a  sharp-nosed 
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body,  like  t  wedge,  is  relatively  far  away  fron  the  surface 
and  the  boundary  layer  at  low  Mach  numbers  except  in  the 
iasiediate  vicinity  of  the  point  (Fig.  1).  However,  at  hyper¬ 
sonic  Mach  numbers,  the  shock  wave  comes  close  to  the  surface 
(Fig. 2).  In  explaining  his  experiments  on  shock  interaction, 
Nagamatsu  comments. 

The  shock  wave  and  boundary  layer  seem  to  be  merged 
before  separating  into  distinct  shock  wave  and 
boundary  layer.  As  the  flow  Mach  number  was  in¬ 
creased,  the  merged  region  extended  over  a  larger 
portion  of  the  p 1  a te ( 1 8 : 461 )  . 

He  also  points  out  some  of  the  variations  in  flow  properties 

this  interactions  canses. 

When  the  shock  wave  and  the  boundary  layer  were 
merged,  the  pressure  from  the  shock  wave  to  the 
surface  for.  a  given  location  from  the  leading  edge 
was  approximately  constant.  After  the  shock  wave 
became  separated  from  the  boundary  layer,  the  pres- 
.  sure  behind  the  shock  wave  was  greater  than  the 
surface  pressure  at  the  same  location  (18:462). 

Boundary-layer  theory  is  not  able  to  predict  any  of  these 
interactive  results.  Therefore,  any  solution  at  the  leading 
edge  of  the  body  probably  has  an  inherently,  large  error. 

This  error  at  the  leading  edge  should  be  ignored  to  get  an 
indication  of  the  overall  error  of  the  boundary-layer  code 
over  the  surface.  The  finite-difference  method  also  has 
problems  predicting  the  large  gradients  at  the  leading  edge. 
Therefore,  neglecting  large  leading  edge  errors  due  to  the 
boundary  layer  solution  also  leads  to  neglecting  a  large 
leading  error  due  to  the  finite-difference  method. 
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PI?.  1  ,  Supersonic  Shock-Boundary  Layer  Interaction 


Pig.  2  Hypersonic  Shock-Boundary  Layer  Interaction 


Chapter  III :  Finite  Difference  Solot ion 
The  bounder y- 1  aye r  equations'*  Eqs  (2)  -(5),  are  a  non¬ 
linear  set  of  partial  differential  equations  which  describe 
flnid  flow  close  to  a  surface.  Their  mathematical  solution 
is  difficult.  However,  using  finite-difference  techniques, 
the  solution  of  the  boundary-layer  equations  and  the  resul¬ 
ting  description  of  the  flow  field  is  possible.  First,  the 
domain  of  the  flow  field  must  be  divided,  or  differenced, into 
a  grid  pattern.  This  grid  may  be  in  any  form  which  makes 
solution  of  the  problem  easiest.  For  this  study,  three  dif¬ 
ferent  grids  which  characterize  three  d  i  f  f  erent,  so  lut  ion 
spaces  are  of  interest:  the  physical  grid,  the  Dorodnitsyn- 
transformed  grid,  and  the  computational  grid.  The  physical 
grid  is  an  orthogonal,  surface-normal  grid  which  describes 
the  flow  pattern  in  the  physical,  compressible  plane.  In  the 
physical  plane,  the  bound ary- 1  ay e r  equations  are  highly  non¬ 
linear  and  vary  with  viscosity  and  density.  Therefore,  the 
Dorodnitsyn  transform  is  introduced  to  take  the  density  de¬ 
pendence  out  of  the  problem.  This  transformation  adjusts  the 
normal  coordinate,  y,  to  put  the  density  dependence  into  the 
definition  of  a  new  normal  coordinate,  T,  instead  of  in  the 
boundary- 1  aye r  equations.  This  creates  coordinates,  X  and  T, 
for  the  Dorodni t syn- t r an s f ormed  grid.  The  final  step  is  to 
take  all  non-linearity  out  of  the  grid  by  transforming  X  and 
T  to  computational  variables,  g  and  q.  Expressing  the  boun¬ 
dary-layer  equations  in  terms  of  the  new,  computational  varia- 
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bles  resalts  in  solation  equations  which  are  aore  easily 
solved  in  a  nnifora,  computational  space.  The  coaputat ional 
(rid  with  coordinates.  ?  and  q ,  has  a  constant,  step  size  of 
one  in  both  the  streaawise  and  normal  directions.  The  physi¬ 
cal  and  Dorodnitsyn  grids  have  varying  step  sizes.  Boundary 
conditions  and  initial  conditions  are  then  added  to  the  solu¬ 
tion  equations,  and  the  equations  are  solved  in  the  computa¬ 
tional  space  using  an  implicit,  Successive-Over-Relaxation 
(SOR)  method. 

Ika  Solution  Grid 

The  solution  grid  divides  the  flow  field  into  exact, 
coordinate  locations.  The  solution  of  the  flow  problem  des¬ 
cribes  the  velocity,  enthalpy,  and  other  flow  conditions 
at  the  grid  locations.  For  each  flow  problem,  there  is  an 
optimum  grid  which  best  resolves  the  flow  field.  The  accur¬ 
acy  of  the  solution  reflects  how  well  the  grid  determines  the 
gradients  in  the  flow  field.  With  very  large  spacing  between 
grid  points,  rapidly  changing  flow  conditions  in  a  particular 
area  may  not  appear  due  to  lack  of  resolution  on  the  grid. 
Increasing  the  number  of  grid  points  increases  resolution  but 
also  increases  computational  time.  It  is  also  a  crude,  brute 
force  approach  to  the  problem.  An  ideal  solution  is  to  put  a 
large  number  of  grid  points  in  the  high-gradient  regions  and 
fewer  points  in  low- gradient  regions  without  increasing  the 
total  number  of  grid  points.  The  boundary- layer  problem 
needs  the  majority  of  the  points  in  the  boundary  layer,  es- 
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pecially  near  the  leading  edge,  where  flow  conditiona  are 
changing  rapidly.  The  least  concentration  of  points  should 
be  in  the  inv isc id—f low  region  ontside  the  boundary  layer 
where  conditions  are  esaentially  constant.  The  particular 
way  this  type  of  grid  is  built  depends  on  the  grid  equation 
chosen. 

This  study  uses  a  one-dimensional,  grid  equation  to  de¬ 
termine  the  grid  spacing  in  one  d imena ion, anal agous  to  a 
multi-dimensional  Poisson  equation  (an  elliptic^  partial  dif- 
ferental  equation).  The  other  dimension  is  then  scaled  using 
a  velocity  profile  of  Blasina',  boundary-layer  solution. 

The  one-dimensional  grid  equations  are 

Xgg+  PXg  -  0  ,  (streamwiae  direction)  or  .  (34) 

*Tpi+  "  0  (normal  direction)  (35) 

where 

{  -  $(x,y, t) 

q  •  q  (  x  ,  y ,  t )  (36) 

P  and  Q  are  control  functions.  The  grid  equations  stretch 
the  grid  points  depending  on  the  control  functions  P  and  Q. 
The  grid  used  to  verify  results  for  Lange's  non- adaptive  grid 
boundary- 1  aye r  program  defines  the  points  in  the  streamwiae 
direction  first.  Setting  P  equal  to  zero  forms  a  constant 
spacing  in  the  x  direction.  Then,  using  a  parabolic  velocity 
profile,  C/0e  ■=  (Y/5)*^  for  6  defined  by  Eq  (8),  determines 
the  y  coordinates.  Due  to  the  constant,  x  spacing  and  the 
parabolic  characteristics  of  the  power  law,  the  power— law 

z2 


grid  is  several,  exp  and  ing,  p  ar  abo  1  ic  curves,  as  Fig.  3  shows. 
The  flow  solution  which  comes  from  using  this  power— law  grid 
can  be  compared  with  the  optimized  grid  solution.  The  opti- 
mized— grid  solution  also  uses  constant  spacing  in  the  x  di¬ 
rection.  However,  it  uses  Eq  (35)  to  fix  the  y  coordinates 
at  an  x  location  and  then  scales  the  rest  of  the  grid  using 
Eq  (8).  Unified  Difference  Relations  (UDR)  given  by  Hodge, 
Leone,  and  McCarty  solve  the  grid  equation  (9).  The  UDR  for 
Eq  ( 3  5 )  a  r e 

Yi+1  ~  <l+e-Q>*i  +  e"Q  Y i„1  =0  (Q  >  0)  (37a) 

eQ  Yi  +  1  -  (eQ+l)Y.  +  Y^j  =  0  (Q  <  0)  (37b) 

A  tridagonal  iteration  scheme  solves  Eqs  (37)  (3:10).  To 
scale  the  grid,  a  ratio  of  6's  at  different  x  locations  is 
set  up  and  like  terms  cancelled.  The  relationship  between 
any  two  x  locations  becomes 

Y2  »  Yj  (X2/Z1)1/2  (38) 

This  study  defines  Y^  and  Xj  and  solves  the  normal  grid  equa¬ 
tion  at  the  far  end  of  the  plate.  For  each  given  X2  loca¬ 
tion,  a  Yj  value  results.  Fig.  4  shows  the  initial  grid  for 
the  optimized  solution.  To  optimize  the  grid,  the  Q  that 
minimizes  the  truncation  error  in  the  third  derivative  of  the 
tangential  velocity  in  the  computational  plane,  U  ^  ,  must  be 
found.  Therefore,  the  grid  must  now  be  transformed  into  the 
computational  plane. 

Although  the  grid  is  computed  in  the  form  presented 
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FI*.  4a  Exponential  Grid 


Fig.  4b  Expanded  Exponential  Grid 


above,  this  grid  is  related  to  coapntational  plane  variables, 
$,  n.  »nd  *•  by  metric  coefficients.  The  functional  rela¬ 
tionships  between  the  transforaed  variables  and  comput a t ions  1 
variables  are 


«  -  $(X.Y,t) 
n  ■  n(X,T, t) 
t  -  t 


(39) 


Using  the  chain  rule  and  solving  for  the  derivative  of  each 
computat ional  variable  with  respect  to  the  transforaed  varia- 


b 1 e  s ,  the 

metric 

coe  fficients  are 

*X  “ 

VJ 

nx  " 

-vj 

*X 

$Y  * 

-vj 

ny  “ 

X?/J 

*Y 

«t  - 

<Vt 

-  XtY,1)/J 

^t  =  (1tT$  ' 

V 

J  - 

X^Yq  - 

V? 

For  this 

study. 

X  does  not 

change 

in  the  normal 

'  \  « 

X^  is  zero.  Also,  for  steady-state  solutions  where 

is  not  changing 

with  time. 

Xt  and 

Y  t  are  zero. 

The 

(40) 


ection,  so 


state  metrics  are 


=  1/X, 


nx  -  -Y5/x?Tr 


nY  =  1/Yn  (41) 


If  the  adaptive  grid  changes  with  time,  the  metrics  also 
include  the  time  dependent  metrics 

5t  =  -xt/x5  nt  =  -Yt/’iY  (42) 

These  metric  relationships.  Eqs  (41)  and  (42),  make  it  pos- 


sible  to  convert  the  boundaxy- layer  equation*  into  computa¬ 
tional  variables.  The  Dorodnitsyn  transformation  changes 
the  physical  grid  shown  in  Fig.  5  to  a  transformed  grid, 
shown  in  Fig.  6.  The  metrics  then  mathematically  change  the 
transformed  grid  into  the  uniform,  rectangular  grid,  shown  in 
Fig.  7.  The  uniform,  computational  grid  makes  solving  the 
boundary-layer  equations  simpler.  It  also  allows  for  finding 
an  optimium  solution  grid  which  is  not  possible  in  the  physi¬ 
cal  plane . 

The  Numer ical .  Differencina  Equations 

Whether  the  equations  are  incompressible  or  Dorodnitsyn- 
transformed,  the  boundary- layer  equations  are  more  easily 
solved  if  they  are  transformed  into  variables  of  the  computa¬ 
tional  plane.  The  constant,  normalized  mesh  size  of  the 
computational  plane  simplifies  the  differencing  equations. 

Eqs  ( 5 )  —  ( 7 )  when  changed  from  X  and  T  in  the  transformed 
plane  to  {  and  q  in  the  computational  plane  become 
Continuity:  Sx®$+1'X^q+,,  =  ®  (43) 

Momentum:  0t  +  U(  +  “  l/Betp^idJ^qylJqy  (44) 

Energy:  Ht  +  D(H^?x  +  Hnqx)  +  ,* 

Pt/p  +  (  JI£_H„i»y>Wy  +  t mlZ£zll  ny(o£Zl),Jn  ny  (45) 

PrRe  .  2RePr 

(15:23).  For  the  momentum  and  energy  equations,  the  differ¬ 
encing  methods  used  are  three-point,  windward  differencing  on 
the  first  derivative  tens,  second-order,  central  differ¬ 
encing  on  second  derivative  viscous  terms  such  as  ,  H_ _ , 

qq  qq 

and  T  ,  and  two-point,  backward  differencing  on  the  time 


terms.  At  the  boundaries,  where  three  points  cannot  approx¬ 
imate  the  first  derivatives,  the  differencing  is  two-point, 
npwind.  This  method  is  first-order  accurate  in  time  and 
second-order  accurate  in  space,  except  at  the  boundaries 
where  it  is  first-order  accurate  in  space.  Lange  uses  this 
method  because,  "this  differencing  scheme  guarantees  diagonal 
dominance  and  is  second-order  accurate  in  the  computational 
plane"  except  at  the  boundary  points  as  noted  above  (15:24). 
After  applying  the  differencing  to  Eqs  (44)  and  (45),  and 
solving  for  U  ^  j  and  H  ^  j  ,  the  solution  equations  result. 
Momentum:  °ni,j  “  +  *5x  AtDuIli-l,j 

-  b?j  AtUDni_2j  +  [cqxO  At  +  cqyV  At 
+  T|y  At((p^i)j_1/2  ny)/Re]Dnit  j_2 
-  tdnxU  At  +  dnYV  At]Dni  j_2 
+  ny  AtC(p#/)j+i/2  ^Y^^i.j  +  l  ,/®o)/DCoef 

where  Ucoef  =  l  +  e?jD  At  +  fiwU  At  +  fnYv  At  + 

HY  At  [Tpji)j+i/2  ^Y  +  j-1/2 

Energy:  11“.,.  =  (Hn-*  i#  j  +  .?xO.  AtH“  j. .  j  -  b?xU  AtH“i.2,j 
+  [cqxD  At  +  cqyV  At  +  nY  A  t  (  (  o/t  )  j  _  1/2  nY)J  ani,j-i 

D  ^  D  .  *  * 


-  (dqxU  At  +  dqyV  pt)Hni>j_2  +  TLY  -AJlH  Pfl  )  j  +  i  /  2  Y^^i.i+l 

Pr  Re 


2  ( 1-Pr ) Re 


C  (  PA4)  j  +  1/2  Y(D“i,j  +  l) 


-  ((p/i)j  +  l/2  *»y  +  <P0>j-l/2  ’)Y)(u“i,j)2 

+  (p//)j_l/2  T»Y(°ni,  j-l>2j)  /  Hcoef  1 

Hcoef  =  1  +  •?XDn'1i.j  At  +  fnx^^i.j  At  +  f^YV"-1!^  At 

+  IL.Y  <  P/1  >  j+1/2  ^Y  +  (p/l>j-l/2  ^Y-l 
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The  i  and  j  give  the  abscissa  and  ordinate  locat ions , respec¬ 
tively,  in  the  compnt a t iona 1  plane.  In  the  equations  above, 
the  coefficients,  a  -  f,  change  depending  on  the  grid  point. 
Using  two-point  backward  or  forward,  or  three-point  forward 
or  backward  differencing  changes  the  coefficients.  Table  I 
has  the  coefficients  for  each  case.  To  linearize  the  eqna¬ 
tions,  Lange  lags  the  the  first  iteration  solution  in  tine  so 
U,  V,  and  H  are  calculated  based  upon  the  U,  V,  and  H  at  the 
sane  location  but  previous  tine  step.  For  subsequent  itera¬ 
tions,  U,  V,  and  H  are  calculated  based  upon  the  U,  V,  and  H 
fron  the  previous  iteration  step.  Next,  the  netric  coeffi¬ 
cients  need  to  be  differenced. 

Fron  the  previous  discussion  of  the  netrics,  the  only 
non-zero  netrics  are  qx,  and  i)y.  The  inverse  relation¬ 

ships  fron  Eq  (41)  generate  the  expressions  for  the  differ¬ 
enced  netrics.  Since  equals  X^  and  a  constant  spacing 
between  streanwise  points  is  used,  a  first-order,  central 
difference  for  Xg  deternines  a  constant  {y.  qy  is  found 

TABLE  I 

COMPUTATIONAL  PLANE  DIFFERENCE  COEFFICIENTS 


i 

j 

a 

b 

C 

d 

e 

f 

2 

2 

1.0 

0.0 

1.0 

0.0 

l.o 

1.0 

2 

3- 

1.0 

0.0 

2.0 

0.5 

1.0 

1.5 

2 

jnax-1 

1.0 

0.0 

-1.0 

0.0 

1.0 

1.0 

3- 

3- 

2.0 

-0.5 

2.0 

-0.5 

1.5 

1.5 

3- 

j  sax-1 

2.0 

-0.5 

-1.0 

0.0 

1.5 

-1.0 

inax-1 

j  nax-1 

2.0 

-0.5 

-1.0 

0.0 

1.5 

-1.0 

similarly  by  using  Y^  ,  which  is  not  constant.  The  spacing  in 
the  normal  direction  does  change  with  each  streamwise  loca¬ 
tion.  Lange  uses  three,  different  schemes  depending  on  the 
location  of  in  the  equation  (15:26).  Y^  in  the  viscoua 

terms  is  differenced  with  a  second-order,  central  scheme 
about  the  j+j/j  or  j-i/2  points.  This  gives  an  overall  sec¬ 
ond-order  accurate  Y^  term.  Y^  for  all  other  terms  is  repre¬ 
sented  by  a  second-order,  central  difference  for  the  interior 
points.  This  study  uses  a  three-point  windward  scheme  at  the 
boundaries.  This  guarantees  second-order  accuracy  for  all  Y^ 
differencing.  For  Yg,  Lange  uses  an  analytic  calculation 
for  the  metric.  Lange  ezplaina. 

Initially.  Yg  was  calculated  using  a  central  dif¬ 
ference  with  a  backward  difference  at  the  end  of 
the  domain.  A  large  amount  of  leading  edge  error 
was  induced  by  this  method  (15:27). 

This  study  uses  an  analytic  metric  different  from  Lange’s. 


However,  since  an  analytic  solution  over  the  entire  plate 
limits  the  method's  general  applicability,  it  is  preferable 
to  use  the  analytic  solution  only  at  the  leading  edge.  Lead¬ 
ing-edge  error  is  inherent  in  the  bound a ry- 1  aye r  solution 
if  Yg  is  numericaly  calulated.  However,  if  Yg  is  calculated 
analytically  with  Eq  (48)  at  the  first  couple  of  streamwise 
locations,  the  numerical  solution  can  then  be  calculated  on 
the  remaining  points  downsteram  of  the  leading  edge  with  less 
error.  Consequently,  Yg  is  calculated  with  Eq  (48)  at  the 
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first  two  streamwise  locations  and  with  numerical,  difference 
equations  at  the  rest  of  the  streamwise  locations.  The  nu¬ 
merical,  difference  equations  use  central  differencing  on 
the  interior  and  three-point,  backward  differencing  at  the 
back  of  the  plate.  This  results  in  Yg  with  the  same  second- 
order  accuracy  as  Y^  and  Xg.  Lange  gives  an  example  of  ti 
resulting  fully-differenced,  linearized,  momentum  and  energy 
equations  (15:28).  Using  Lange's  SOR  method  to  compute  U  and 
H  from  the  differenced  equations  gives  values  used  in  the 
continuity  solution  for  V  (15:30-32). 

The  transformed,  continuity  equation,  Eq  (43),  now  con¬ 
tains  only  one  unknown.  V  has  not  been  solved.  Expressing 
Vn  in  terms  of  the  known  metrics  and  U,  V  becomes 


+  h-\ 

S 


(49) 


Lange  integrates  Eq  (49)  with  the  trapezoidal  rule  for  a 
constant  step  of  one.  Finite  differencing  on  Ug  is  a  three- 
point,  windward  scheme  (15:32).  The  integration  of  U^  used 
in  this  thesis  is  different  than  Lange's  treatment.  The 
term  separates  into  two  parts. 

Y*  U  =  1_[Y£  U  -  fu(Y«)  dq  ]  *  (50) 

*5  .  j-1 


The  integral  is  evaluated  using  a  trapezoidal  rule  with  Y^ 
averaged  about  the  j-1/2  point  prior  to  the  integration. 
After  collecting  terms,  the  difference  equation  for  the  in¬ 
terior  points  is 


34 


‘I\ 


Vi.j  *  Vi,j-1  *  t<Y5i.j  + 

+  <Yi.j  '  Yi. j-l>*C-°*75<Ui. j  +  <51> 

+  1*0<Di-l.j  +  °i-l,j-l>  “  °-2*<*i-2,j  +  Wi-2.  J-1>1  J 

Solving  for  V  completes  the  flow  solution.  All  parameters  of 
the  boundary  layer  problem  are  then  known  or  can  be  deter¬ 
mined  from  U,  V,  and  H. 

To  determine  how  close  the  computed  boundary  layer  solu¬ 
tion  duplicates  the  heat  transfer  coefficient,  h,  Eqs  (26) 
and  (27)  are  solved  numerically.  Since 


k '  *  fi’C  /Pr 


the  transformed,  heat  equation  becomes 


(52) 


-  U'C  'o  '  /dT'\ 

Pr  P  *  (dY '  )y=0 


(53) 


The  dT’/dY'Jy^Q  term  is  non- d imens iona 1 i z ed  and  Eq  (53)  is 
then  transformed  to  the  computational  space  (15:33).  ,  This 
equation  is 


PrLRe  U  q7 


(54) 


where  dT/dq  =  (-3  T.j  +  4  Ti>2  ~  Ti3)/2  and 

qY  =  l/Y,,  -  2/(-3Yi>1  +  4Yi>2  -  Yl>3) 

The  equation  differencing  is  three-point  upwind  for  the  me¬ 
tric  Yjj  and  dT/dq,  Eq  (26)  is  rearranged  so  that 


h'  =  q  '  /  (T  ’ -T  ') 

*i»  ' 


(55) 
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(22:95).  The  heat  transfer  coefficient  calculated  ii  then 
compared  to  the  theoretical  h  from  Eq  (32).  If  the  numerical 
scheme  is  accurate  the  ratio  of  the  two  h's  should  be  1.0  for 
incompressible  cases. 

Another  parameter  that  is  calculated  to  compare  the 
computer  solution  to  theoretical  incompressible  results  is 
Cf.  Eq  (9)  gives  the  Blasius  solution  for  the  skin  friction 
coefficient.  This  relation  is  derived  from  the  definition  of 
Cf  which  is 


The  non-d imens iona 1  i z e d ,  transformed  equation  is 


Cr  = 


(57) 


Y=0 


Three-point,  backward  differencing  evaluates  and  Y^  at  the 
wall.  When  the  computed  from  Eq  (57)  is  divided  by  the 
theoretical  C £  from  Eq  (9),  the  best  solution  produces  a 
ratio  closest  to  1.0  for  the  incompressible  case. 

Boundary  Cond i t ions 

The  flow  conditions  at  the  wall  and  the  upper  boundary 
of  the  domain  fix  the  boundary  conditions.  Assuming  no  slip 
conditions  at  the  wall,  uw  and  Uw  are  zero.  Also,  there  is 
no  flow  through  the  boundary  surface,  therefore,  vw  and  Vw 
are  zero.  The  surface  temperature,  Tw ,  for  the  isothermal 
wall  cases  is  530°R,  Freestream  temperature  is  also  530  °R 
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for  most  of  the  incompressible  runs.  These  cases  come  close 
to  an  adiabatic,  wall  condition.  At  the  domain's  npper  boun¬ 
dary,  the  edge  conditions  depend  on  the  given  Mach  number. 

For  subsonic  flows,  the  edge  conditions  are  the  same  as  the 
freestream  conditions.  At  Mach  1  and  above,  oblique  shock 
wave  theory  sets  the  flow  conditions  behind  the  shock.  The 
conditions  behind  the  shock  are  the  edge  conditions  for  su¬ 
personic  flows.  Unlike  these  fixed  boundary  conditions,  the 
initial  conditions  change. 


Initial 


For  the  first  iteration  of  the  numerical  solution,  an 


initial  approximation  of  U,  H,  and  V  must  be  assumed.  A 


cubic,  velocity  profile  is  a  relatively  simple  approximation 
to  the  boundary  layer  and,  as  shown  in  Schlichting,  it  comes 
close  to  Blasius',  exact,  boundary- 1 ay e r ,  velocity  profile 
(23:206).  The  cubic,  velocity  profile  is 


u  =  u0  *  [1.5  (  y  *  /  5  )  -  .5  (y'/S)3]/ua, 


(58) 


For  compressible  solutions,  Y'  replaces  y’  and  A(X)  replaces 
6(x)  in  Eq  (58).  The  enthalpy  profile  is  derived  from  a 
temperature  profile  which  is  an  approximate  solution  of  the 
energy  equation  (22:34). 


i.5  ixii  -  . 5i7 v ’ a3 
(8t)  If  8  t ) 


(59) 


Since  the  total  enthalpy,  H,  is  CpT+U3/2  for  isentropic 
flows,  Eq  (55)  yields  the  non-d ime n s  i  ona 1 i z ed  gue s s  for  H. 
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Hi.j  -  Hw  + 


(He-Hw) 


(60) 


The  thermal,  bounds ry- 1  aye r  thickness,  &t,  in  Eq3  (59)  and 
(60)  is 


6t  -  6/(1.026  Pr1/3) 


(61) 


This  value  for  8t  is  valid  if  the  entire  plate  surface  is 
heated  and  if  6t  is  zero  at  the  front  of  the  plate  (22:35). 
For  Pr*l,  8t  will  be  less  than  8.  Finally,  initial  v  values 
are  assumed  to  be  zero  throughout  the  domain.  These  complete 
the  initial  guess  for  the  flow  solution.  The  guess  is  input 
along  with  the  grid,  and  solution  of  the  difference  equations 
gives  a  final,  numerical  solution  for  the  flow  characteris¬ 
tics.  Unfortunately,  the  input  grid  may  not  be  able  to  ade¬ 
quately  resolve  some  of  the  gradients  in  the  solution. 
Therefore,  an  optimized  grid  which  minimizes  the  truncation 
error  and  produces  better  flow  solutions  must  be  found. 


Using  tinitb-difference  methods  to  solve  the  boundary- 
laver  equations  requires  polynomial  approximation  of  the 
first  derivatives  of  velocity  and  enthalpy.  The  polynomial 
approximations  are  not  exact.  Truncation  error  is  inherent 
in  the  solution  since  all  terms  of  the  approximation  are  not 
included  in  the  calculations.  In  the  case  of  second-order 
differencing  of  velocity  derivatives,  truncation  error  is 

some  multiple  of  the  third  derivative  terms  such  as  U _  and 

Um.  Minimizing  the  third  derivative  terms  insures  the 
truncation  error  is  small,  and  the  velocity  solution  is  con¬ 
sequently  most  accurate.  Since  the  velocity  gradients  in 
the  streamwise  direction  are  small  compared  to  the  velocity 
gradients  in  the  normal  direction,  this  study  minimizes  only 
U^  terms.  Also,  as  U^^  approaches  zero,  approaches 

a  constant  and  U^  approaches  zero.  Finding  the  solution 
grid  which  minimizes  U  ^  thus  optimizes  the  boundary- layer 
solution.  Powell's,  conjugate-direction  method  which  iter¬ 
ates  on  an  n-d imens iona 1  array  until  finding  the  minimum  of  a 
desired  function  without  calculating  the  derivatives  of  the 
function,  is  well-suited  to  this  problem.  Powell’s  method 
finds  the  grid  control  function,  Q,  in  Eq  (35)  which  forms 
the  optimum  grid. 

Truncation  error  is  present  in  any  finite-difference 
solution  which  uses  polynomial  approximations  to  define  gra- 
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dient  terms.  The  finite  differencing  of  the  boundary- layer 
equations,  uses  Taylor,  polynomial  expansions  of  the  gradient 
terms.  For  example,  the  first-order,  three-point,  backward 
difference  of  a  gradient  term,  is 


i.j 


-4.  J! 


i , 


W,  4-2  +  1  +  1,  +...  (62) 

—  i.J-2  -  qqq  - 


assuming  Aq  is  1.  Since  the  first  three  terms  replace  a 

gradient  term,  like  ,  in  the  f in i t e- d i f f e r e nc ed ,  boundary- 

layer  equations,  the  truncation  error  is  the  sum  of  all  the 

higher  derivative  terms.  The  most  heavily  weighted  term  of 

the  truncation  error  is  the  term.  If  W___  is  zero, 

W__„_  is  zero.  Therefore,  minimizing  the  third  derivative 
1111 

terms  should  significantly  decrease  the  truncation  error  in 
the  finite-difference  solution. 

The  boundary- layer  equations  are  two-dimensional  equa¬ 
tions.  For  a  completely,  accurate  solution,  two-dimensional, 
grid  equations  should  define  the  solution  grid.  However, 
this  study  uses  the  characteristics  of  the  incompressible, 
boundary  layer  to  simplify  the  grid  selection.  One  of  the 
bounda ry- 1  aye r  simplifications  of  the  Navier-Stokes  equations 
is  that  the  normal  velocity  gradients  are  much  greater  than 
the  streamwise,  velocity  grad'euts.  Therefore,  any  error  in 
the  normal  velocity  gradients  has  a  larger  impact  on  the 
overall  error  of  the  computer  solution  than  error  in  the 
streamwise  gradients.  Making  the  x  locations  constant  in 
the  adaptive  grid  allows  the  grid  to  concentrate  on  optimiz- 
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ing  only  the  gradients  in  the  normal  direction.  Primarily 


* 

the  adaptive  grid  tries  to  minimize  the  sum  of  U^'  ,  bat  it 

can  also  minimize  the  sums  of  D„„  ,  U_  ,  or  some  other  com- 
bination  of  these  terms.  Once  the  normal  direction  is  opti¬ 
mized,  the  ttreamvise  dimension  is  scaled  according  to  Bias- 
ins'  description  of  the  incompressible,  boundary  layer.  A 
parabolized,  grid  scheme  can  be  nsed  to  avoid  this  scaling. 

Eq  (38)  shows  a  sqnare-root  dependence  on  each  streaawise 
coordinate  of  the  Blasins  solution.  Ordinate  locations  found 
using  a  given  grid  control  function  and  the  normal  grid  equa¬ 
tion,  Eq  (33),  at  one  x  location  are  scaled  bn  the  rest  of 
the  grid  using  Eq  (38).  Powell's  minimization  method  finds 
the  optimum,  grid  control  function,  Q,  which  produces  an 
optimized  flow  solution. 

Powe 1 1 ' s  Method 

Powell's,  minimization  method  has  several  features  which 
make  it  attractive  for  optimizing  the  normal  gradients. 
Powell's  method  does  not  require  calculating  derivatives  of 
the  function  being  minimized.  This  is  especially  helpful  in 
this  study  where  no  explicit  function*  is  being  minimized. 
Powell's  method  is  flexible.  Any  computed  quantity  can  be 
minimized.  Any  parameter  can  become  the  minimized  function 
in  Powell's  algorithm.  This  allows  the  optimization  of  the 
grid  with  several,  different,  minimized  parameters  to  find 
the  best  optimization  technique,  Powell's  method  assumes  a 
locally,  quadratic  form  between  points  rather  than  a  linear 


fora.  Powell's  method  also  converges  quickly  i&  a  finite 

number  of  steps.  Powell  explains, 

when  the  procedure  is  applied  to  a  quadratic  form 
...the  ultimate  rate  of  convergence  is  fast  when 
the  method  is  used  to  minimize  a  general  function 
(19:155). 

Powell's  use  of  two  conjugate  directions  to  search  for  the 

minimum  achieves  this  fast  convergence.  Therefore,  Powell's 

method  is  accurate,  allows  fast  convergence,  is  flexible,  and 

does  not  require  further  derivatives  of  U„ ....  to  minimize  this 

nnn 

parameter. 

This  study  adapts  Powell's  method  to  minimize  the  sum  Of 

7.  2  7 

0nnn  ,  Unn  ,  0^  ,  or  any  combination  of  these  sums  over  the 

entire  solution  grid.  One  iteration  of  Powell's  minimization 

procecure  is  explained  below. 

(i)  For  rsl,2,...,n  calculate  Xr  so  that 
f(Ar_j  +  XrEf)  is  a  minimum  and  define 

Ar  =  Ar-1  +  ^r®r  * 

(  i)  Find  the  integer  m,  1  <.  m  <  n,  so  that 
( f ( A  j ) - f ( Am) } i s  a  maximum,  and  define 
A=f(Xa_1)-f(AB). 

(iii)  Calculate  fn=  f(2An~A0),  and  define 
fx“  f ( Aq )  and  f2=  f(An). 

(iv)  If  either  f32.fi  and/or 

(fx  -  2  f  2  +f3)*(f1-f2~A)2  >  .SAtfj-fj)2  nse  the 
old  directions  Ej  ,  Ej , . . . ,  En  for  the  next  itera¬ 
tion  end  use  AQ  for  the  next  Aq,  otherwise 

(v)  defining  E=  (An~A),  calculate  X  so  that 

f(A&  +  XE)  is  a  minimum,  use  Ej ,  E2 , . . . , E  2 ,  Effl, 

Em+j,  Effl+2 , . . . En , E  as  the  directions  and  A  +XE  as 
the  starting  point  for  the  next  iteration  (19:156). 

For  this  study,  the  sum  of  O^^2  ,  U^2,  U^2,  or  *  combina¬ 
tion  of  any  of  these  terms  is  the  above  function,  f.  The  Aq 
to  A  parameters  above  are  A(,),  A(2),...,  A(n).  Q,  the  con- 


trol  function,  is  some  function  of  the  A(l)  to  A(n)  parame¬ 
ters. 

Q  -  Q  (A(l) ,  A(2) . . . . ,A(n) )  (63) 

The  X  directions  that  Powell  mentions  are  the  directions 
the  algorithm  searches  to  find  the  function's  minimum.  Ini¬ 
tially,  this  study  uses  two  parameters  so  that 

Q  -  A ( 1 )  +  A( 2  )  (64) 

The  first  step  in  the  minimization  procedure  is  finding  f(Ag) 
which  is  also  f^.  The  flow  chart  in  Appendix  C  shows  how 
program  MNTRER  in  Appendix  D  applies  Powell's  algorithm. 

The  program  MNTRER  in  Appendix  D  finds  each  Xr  by  iter¬ 
ating  over  k  iterations  until  it  reaches  a  minimum  of 
f(Ar_j  +  XrEf) .  For  each  r,  A  is  initially  zero.  Two  other 
values  of  X,  A+  and  X-,  are  also  chosen  a  fixed  percentage 
distance  from  X.  The  k  iteration  then  begins.  Three  new 
parameters,  A+,  A",  and  Ab ,  are  calculated  for  each  r  with 
X  ,  X  ,  and  X . 

Ar<  >  -  Ar(  >  +  Ar.Er  (65) 

The  new  A  values  add  according  to  Eq  (64)  and  form  Q+,  Q~ ,  , 
and  .  Three  new  grids  axe  then  generated  tor  the  + ,  - ,  and 
b  cases.  Three  sets  of  bound a ry- 1  ay e r  equations  are  solved. 
Taking  the  three  b ound a ry- 1  ay e r  solutions,  subroutine  NEWP1 
calculates  0^  ^  an^  *n^  three  functional  values  of 

f,  f+,  and  f~ .  Subroutine  NEW"  also  fits  a  quadratic  to 


X+,  X- ,  and  X.  It  then  finds  the  minimum  of  the  quadratic  to 
produce  a  new  X.  Eq  (65)  then  finds  a  new  A  .  From  this  a 
new  Q  is  found,  and  the  solution  grid  is  generated.  Then, 
the  boundary- 1  aye r  equations  are  solved,  and  fj^  is  deter¬ 
mined.  The  minimum  of  f+ ,  f ~ ,  f,  and  fN  is  determined.  The 
X  associated  with  the  minimum  becomes  the  new  X  for  the  next 
k  iteration.  A  new  value  of  Ar  is  also  found  with  the  new 
X.  Then,  the  k  iteration  increments  and  continues  until 
reaching  preset  maximum  k  or  the  difference  between  the  new 
and  old  X  is  sufficiently  small. 

After  the  k  iteration,  Powell's  step  (i)  is  cc;.lete. 
The  program  then  continues  with  the  other  steps.  Step  (ii) 
compares  the  f  values  for  each  final  Ar  and  determines  the  r 
which  maximizes  f(Ar_^  -Ar).  Program  MNTRER  then  finds  fj 
in  Powell's  step  (iii).  The  program  then  performs  the  tests 
in  step  (iv).  If  the  tests  are  true,  the  Er  values  of  the 
first  iteration  remain  the  same  and  another  iteration  is 
run.  If  the  tests  are  false,  step  (v)  calculates  a  new  value 
of  E  equal  to  An-Ag.  Using  the  same  steps  as  the  previous  k 
iteration,  the  program  finds  the  X  which  minimizes  f(An+XE). 
An+XE  then  replaces  Aq  as  the  starting  point  for  the  next 
fnll  iteration  of  the  procedure.  The  new  E  calculated  in 
step  (v)  also  replaces  the  Er  in  which  the  largest  decrease 
was  made  in  the  last  iteration.  The  whole  procedure  is  re¬ 
peated  until  the  changes  between  Af's  are  sufficiently  small 
or  the  maximum  number  of  iterations  is  reached.  Powell's 
optimization  is  then  complete. 


Powell's  optimization  method,  as  adapted  in  program 
MNTRER,  provides  some  flexibility  to  test  which  minimized 
function  gives  the  most  accurate  results.  Since  Powell's 
method  minimizes  the  function,  f,  any  function  f  can  be  cho¬ 
sen.  This  study  checks  the  affect  of  minimizing  other  deriv¬ 
atives,  or  a  combination  of  derivatives,  on  the  accuracy 
of  the  final,  optimize-i.  boundary-layer  solution.  The  sum  of 
0^  ,  ,  or  any  combination  of  these  three  can  be 

the  minimized  function,  f.  These  are  the  three  functions 
this  study  investigates,  but  any,  definable  parameter  can  be 
minimized  by  the  optimization  method.  Based  on  the  trunca¬ 
tion  error  argument,  minimizing  with  Powell's  method 

should  produce  the  most,  accurate,  f in i t e-d i f f e r enced ,  boun¬ 
dary-layer  solution. 


S 
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Chant  e  r  V :  Results  and  Discussion 


This  study  shows  how  an  optimized,  adaptive-grid  solu¬ 
tion  of  a  boundary-layer  code  produces  a  better  description 
of  flow  characteristics  in  the  boundary  layer.  Before  look¬ 
ing  at  the  optimization  of  the  boundary-layer  code,  the 
characteristics,  capabilities,  and  limitations  of  the  boun¬ 
dary-layer  code  must  be  understood.  For  example,  questions 
which  must  be  answered  are  how  do  initial  conditions  affect 
the  solution,  or  how  many  iteration  and  time  steps  are  neces¬ 
sary  to  get  a  converged  solution.  The  boundary-layer  program 
used  to  solve  the  f in i t e-d i f f e r e nc ed .  boundary- layer  equa¬ 
tions  has  several  improvement's  to  Lange's  boundary-layer 
code.  Comparison  of  the  computed,  boundary- layer  solution 
and  Blasius',  theoretical  solution  for  incompressible  flow 
over  a  flat  plate  shows  the  increased  accuracy  of  the  new, 
bounda r y- 1  ay e r  code.'  Powell's  optimization  method  is  then 
applied  to  the  new,  bound ary- 1  ay e r  code  to  optimize  the  adap¬ 
tive  grid.  The  optimized,  adaptive  grid  is  more  accurate 
than  the  non-adaptive  grid,  incompressible  solution.  The 
effect  of  different  parameters  on  the  optimization  method  is 
investigated.  These  investigations  of  the  incompressible 
cases  verify  the  new,  bound a r y- 1  ay e r  code  and  the  optimiza¬ 
tion  method.  Then,  the  optimized,  adaptive  grid  is  applied 
to  supersonic  and  hypersonic,  compressible,  flow  problems. 
Solutions  for  compressible  flow  over  the  flat  plate  and  wedge 


•re  compared  to  previous  theoretical  development.  This 
comparison  shows  the  applicability  of  the  optimized,  adaptive 
grid  and  the  boundary-layer  code  to  compressible,  flow  prob¬ 
lems.  The  first  step  however  is  to  verify  that  the  bonndary- 
layer  code  can  solve  for  the  correct,  incompressible  resnlts. 

A3?, Iff  Plow 

The  first  step  in  determining  the  applicability  of  the 
optimisation  method  verifies,  that  the  code  reproduces  the 
exact,  theoretical  results  for  incompressible  flow  over  a 
flat  plate.  The  boundary-layer  oode  is  set  up  to  handle 
compressible,  flow  problems  where  density  and  viscosity  are 
changing.  The  input  flow  conditions  need  to  model  incompres¬ 
sible  flow  where  density  and  viscosity  should  be  essentially 
constant.  To  do  this,  the  incompressible  cases  are  run  at 
Mach  .01.  Freestresm  temperature  and  pressure  are  set  at 
530°R  and  2116.2  psf.  Surface  temperature  is  also  a  constant 
530°R.  For  these  conditions  density  and  viscosity  are  essen¬ 
tially  constant.  Throughout  the  runs,  the  non-dimensional- 
ized  density  ranges  from  .99941217  to  .99941725  in  the  nor¬ 
mal,  flow  direction.  The  density  computed  at  the  edge  boun¬ 
dary  is  .99941725.  The  non-d imeus ional ized  density  should  be 
1.0  for  the  entire  domain  and  especially  at  the  edge  boun¬ 
dary.  The  computed  p  is  5,8  X  10”^  %  in  error.  This  is 
insignificant,  but  the  error  is  present.  Lack  of  more  sig¬ 
nificant  figures  in  the  gas  constant  used  and  round-off 
error  in  the  computer  cause  this  error.  The  computed  den- 
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sity  changes  less  than  5  x  10~^%  over  the  domain.  The  densi¬ 
ty  is  therefore  essentially  constant.  The  viscosity  ranges 
from  3.7989885  x  10-7  to  3.7989486  x  10“7,  the  freestream 
valne.  With  this  error  of  1  x  10“^%,  viscosity  is  also 
essentially  constant.  The  inpnt  freeatream  conditions  cause 
the  boundary- layer  code  to  calculate  essentially  an  incom¬ 
pressible  problem.  To  compare  various  characteristics  of  the 
boundary- 1  aye r  code,  a  61x30  solntion  grid  is  used,  except  as 
noted. 

Boundary  Laver  Code 

Convergence  Parameters 

Several,  input  parameters  control  whether  the  boundary- 
layer  solution  converges  on  a  final  solution  and  how  fast  it 
converges  on  that  solution.  These  parameters  are  the  number 
of  time  iterations  (NT),  the  number  of  iterations  at  each 
time  step  (KT) ,  the  size  of  the  time  step  (DT),  the  conver¬ 
gence  criteria  (EPS),  and  the  initial  conditions  for  the 
solution.  Subroutine  BLIMP,  listed  in  Appendix  D,  represents 
the  variables  as  NT,KT,DT,  and  EPS,  respectively.  For  a 
particular  time  step,  the  solution  is  converged  when  the 
maximum  of  the  differences  of  D,  V,  or  H  at  the  old,  itera¬ 
tion  level  and  the  new,  iteration  level  is  below  a  specified 
error  value.  Picking  a  smaller  minimum  error  tolerance,  EPS, 
results  in  a  more  precise  answer.  However,  it  also  means  the 
solution  will  need  more  iterations  and  computer  time  to  con¬ 
verge.  For  this  study,  EPS  is  1.0  x  10~®.  This  gives  ac¬ 
ceptable  precision  without  unnecessary  use  of  computer  time., 
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The  ratio,  Cf/Cf^,  listed  in  Table  IIB  shows  how  close  the 
competed  solution  comes  to  Blasius’,  theoretical  solution. 

TABLE  lib 

Non-Adaptive  Incompressible  Grid  Solutions-  Results 


Cf/Cf t 

!  ad  ins  < 


.91378 

.97052 

.98740 

.99203 


HIM] 


.98827 
.99332 
.99370 
99373 


.97052 
.98740 
.  99203 
99328 


.91378 
99328 


.97051 

.98740 

.99203 

,99328 


Cf/Cf. 


BELT 
-4 


DELE' 
-4 


90550 

4,840 

4.1 

97400 

5.5 

.117 

98961 

H 

VO 

.185 

99193 

.916 

.0277 

99229 

.168 

.00592 

99231 

7.84 

.00996 

99239 

.742 

.00625 

99241 

.0505 

.00433 

99241 

.00199 

97398 

55.3 

1.17 

98971 

6.1 

.186 

.99193 
22 


.90550 

99228 


.97374 

.98956 

.99193 

.99229 


.  918 

,168  _ 


342.0 

.119 


4821.8 

55.67 

6.18 

.1685 


.027 

,00 


8.99 

,000688 


40.9 

1.178 

.1887 

.00060 


*  :  DELE=  maximum  difference  in  U,  V,  or  H  between  itera¬ 
tion  steps:  DELT=  maximum  difference  in  U,  V,  or  H  between 
time  steps 
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The  numbers  shown  are  values  of  Cf/Cft  at  the  x  locations  one 
point  inside  the  leading  edge  and  at  the  trailing  edge. 

Figs.  8-12  show  the  trend  for  how  Cf/Cf{  changes  in  the 
streamwise  direction.  Flow  properties  at  the  first  two 
points  inside  the  leading  edge  are  calculated  analytically 
and  are  not  changed  by  the  numerical  solution.  These  analyt¬ 
ical  points  are  theoretically  the  closest  the  solution  can 
come  to  the  flow  properties.  Leading  edge  error  distorts 
the  results  after  the  analytically  calculated  points.  How¬ 
ever.  as  the  solution  varies  in  the  streamwise  direction, 
the  numerically  calculated  solution  recovers  to  its  true 
value  toward  the  back  of  the  plate.  Therefore,  the  first 
point  shows  the  analytic  value  of  the  ratios,  and  the  trail¬ 
ing,  edge  value  shows  the  best,  numerical  solution  for  the 
run.  Two  other  results  shown  are  DELT,  the  maximum  differ¬ 
ence  of  U.V,  or  H  between  time  steps,  and  DELK,  the  maximum 
difference  for  0,  V,  or  H  between  iteration  steps.  DELT 
shows  the  improvement  in  the  solution  between  successive  time 
steps.  If  the  maximum  difference  for  U,  V,  or,  H  between 
iterations  steps,  DELK,  is  below  1x10“**,  the  solution  is 
converged.  The  comparison  of  values  in  Table  II  shows  the 
effect  of  the  input  parameters  on  the  bound a ry- 1  ay e r  solu¬ 
tion. 

A  parameter  chosen  to  reduce  computer  time  is  the  size 
of  the  time  step,  DT.  For  each  iteration  of  time,  the  time 
in  the  solution  increases  by  DT.  Since  the  boundary-layer 
equations  are  solved  implicitly  by  SOR,  large,  time  steps  do  1 


51 


not  cause  any  instability  in  tbe  solution  of  the  boundary- 
layer  equations.  The  problem  will  still  converge  on  to  the 
steady-state  answer,  regardless  of  the  size  of  the  tine 
step.  Therefore,  to  get  a  converged  solution  in  as  few 
time  steps  as  possible,  this  study  uses  a  large,  non-dimen¬ 
sional  ized , t ime  step  of  1.0x10®.  In  runs  16-20,  the  solu¬ 
tion  still  converges  in  20  time  steps  using  DT  equal  to  100. 
Larger  DT  produces  better  convergence,  however.  DELK  and 
DELT  go  down  as  higher  DT's  are  used.  DELT  and  DELE  in  runs 
1-5  is  smaller  than  the  same  two  parameters  in  runs  10-13, 
and  runs  10-13  are  smaller  than  runs  16-19.  Other  than  this 
improvement  in  DELT  and  DELE,  increasing  DT  does  not  change 
the  solution  much.  Therefore,  to  be  guaranteed  the  best, 
converged  solution  possible  the  larger  DT  is  used  in  the 
steady- state  solutions. 

The  number  of  iterat ions  per  time  step,  ET,  determines 
how  accurate  the  results  are  at  each  time  step.  If  the 
number  of  iterations  per  time  step  is  large,  i.e.  10  in  runs 
6-9,  the  program  does  more  iterations  of  the  equations  per 
time  step.  The  question  arises  whether  it  is  better  to  have 
more  iterations  per  time  step  and  fewer  time  steps  for  con¬ 
vergence  or  fewer  iterations  and  more  time  steps.  Runs  6-9 
set  ET  at  10.  Runs  1-5  set  ET  at  4.  As  expected,  conver¬ 
gence  at  each  time  level  listed  is  better  than  for  corres¬ 
ponding  time  steps  of  runs  1-5.  With  ET  at  10,  better 

results  for  Cf/Cft  are  obtained  at  15  time  steps  (run  8)  than 
are  obtained  after  20  time  steps  with  ET  of  4  (run  5). 


However,  IS  time  steps  with  ET  at  10  takes  37, 5  CPU  seconds 
on  the  CYBER  £QZQ •  For  ET  at  4,  20  time  steps  takes  30  CPU 
seconds.  This  indicates  there  are  more  iterations  of  the 
solution  with  the  lairger  ET.  More  iterations  produce  better 
results  for  Cf/Cft,  DELE,  and  DELT.  Since  a  large  DT  ensures 
the  solution  is  converged,  the  total  number  of  iterations  is 
the  important  factor.  Whether  the  iterations  are  mostly  with 
each  time  step  or  with  more  time  steps,  the  converged  solu¬ 
tion  occurs  when  enough  total  iterations  have  been  made. 
Therefore,  since  acceptable  convergence  is  achieved  with 
/  less  computer  tipe  for  smaller  ET  values  in  runs  1-5,  the 

inputs  for  run  5  are  used  as  a  set  of  baseline  inputs  for  the 
rest  of  the  study. 

Program  Logic  Choices. 

In  this  study,  three  major  changes  have  been  made  to 
Lange's,  boundary- 1  aye r  code.  The  initial  conditions  are 
slightly  different,  the  integration  of  the  continuity  equa¬ 
tion  is  improved,  and  some  of  the  metrics  are  solved  more 
precisely.  These  changes  help  the  boundary-layer  program  to 
more  closely  conform  to  boundary- layer  theory.  The  improve¬ 
ment  is  seen  in  each  case  by  comparing  the  ratio  of  the 
computed  and  theoretical  values  for  Cj.  The  closer  this 
ratio  is  to  1.0  indicates  better  accuracy.  Therefore,  each 
of  the  changes  is  compared  to  find  the  best  method  to  get 
closest  to  Blasius’,  exact,  incompressible,  boundary— 1 ayer 
so lut ion . 

The  choice  of  initial  condit  ions  depends  on  the  physics 
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of  the  flow  problem.  However,  the  f in  1 1 e- d 1 M e r e nc e  solution 
converges  on  the  same  answer  regardless  of  the  initial  condi¬ 
tions  chosen.  The  boundary- layer  code  changes  Lange's,  ini¬ 
tial  guess  for  the  enthalpy  profile,  the  thermal,  boundary- 
layer  thickness,  and  the  U  and  H  values  given  to  the  jmax+1 
points.  Although  the  points  outside  the  jmax  location  are 
not  of  interest  to  the  solution,  the  finite-difference  scheae 
requires  that  points  at  jmax+1  be  defined.  In  Lange's  pro¬ 
gram,  these  jmax  +  1  values  of  IJ  and  H  are  set  at  0.0.  Since 
these  points  are  on  the  boundary,  well  outside  the  boundary 
layer,  the  jmax+1  points  should  equal  the  jmax  values.  These 
values  are  equal  to  the  edge  conditions.  Therefore,  At  jmax 
and  jmax+1  positions,  0  and  H  are  defined  as  0„  and  J. . 

Lange  also  assumes  the  thermal,  b ound 8 ry- 1  ay e r  thickness,  6 t , 
equals  the  bound a ry- 1  ay e r  thickness,  6.  This  approximation 
is  acceptable,  but  it  is  more  accurate  in  the  hypersonic 
limit  and  for  flows  with  Pr  equal  to  1.  However,  for  the 
incompressible  case,  Eq  (59)  is  accurate  in  the  more,  general 
case  where  Pr  is  non- unity  and  the  surface  is  uniformly 
heated  along  its  entire  length.  This  study  .. «  *  s  this  equa¬ 
tion  to  also  form  Eq  (58),  a  better  guess  of  the  enthalpy 
profile.  Lange's  guess  does  not  include  the  last  two  terms 
of  Eq  (58).  These  three  changes  to  the  initial  conditions  do 
not  affect  the  final  solution.  Bo und a ry- 1  ay e r  solutions  with 
all  other  inputs  constant  but  different,  initial  conditions 
converge  to  the  same  final  solution.  Runs  1-5  and  14-15 
compare  the  results  of  the  two  runs.  The  new,  initial  condi- 
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tions  decrease  tlie  accuracy  of  the  solution  at  the  first  tine 
step.  But  the  final  solution  is  slightly  more  accurate  with 
the  new  initial  conditions.  Although  it  is  not  necessary  to 
know  the  physics  of  the  problem  fully  to  get  the  correct 
initial  conditions,  accurate  initial  conditions  increase  the 
accuracy  of  the  boundary- layer  solution  by  decreasing  the 
iterations  need  for  convergence. 

Another  change  to  Lange's  program  is  in  the  integration 
of  the  continuity  equation  to  solve  for  V.  Eq  (51)  shows  the 
integration  of  .  Fig.  8  shows  the  dramatic  increase  in  the 
accuracy  of  the  new,  boundary- 1  aye r  code.  Both  runs  are 
converged.  The  solution  shown  using  Lange's  integration 
method  for  uses  a  two-point,  windward  difference  at  the 
boundaries  with  the  two  leading-edge  points  defined  analyti¬ 
cally.  The  new,  boundary- layer  code  uses  three-point,  wind¬ 
ward  difference  with  three  leading  edge  points  defined  ana¬ 
lytically.  This  accounts  for  the  large  overshoot  for  Lange's 
solution,  while  the  new  code  undershoots  slightly.  The 
important  result,  seen  in  Fig.  8,  is  the  new  code  recovers  to 
more,  accurate  values  of  Cf/Cf^.  Lange's  code  approaches 
.93319,  while  the  new  code  approaches  .99229.  The  errors  are 
6.7%  and  .71%,  respectively.  The  order  of  magnitude  increase 
in  accuracy  for  Cf/Cf^  indicates  he  new,  bounda ry- 1  ay e r  code 
is  better  at  predicting  the  incompressible  solution. 

As  mentioned  earlier,  the  new,  boundary- layer  code  uses 
a  different,  finite-difference  appr os ima t ion  of  the  metrics, 
and  1^.  Using  the  three-point  windward  scheme  at  the 
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boundaries  gives  second-order  accuracy  throughout  the  netrics 
solution.  One  differencing  method  tried  uses  two-point, 
windward  differencing  which  is  only  first-order  accurate. 

For  the  Xg  metric,  the  differencing  does  not  change  the 
solution.  Since  the  streasiwise  grid  spacing  is  constant,  , 
either  differencing  scheme  gives  the  same  values  of  Xg.  The 
differencing  on  Yg  does  matter,  though.  Using  the  new'  inte¬ 
gration  of  ,  solutions  with  three-point  and  two-point, 
windward  differencing  on  Yg  are  run.  Fig.  9  shows  the  ratio 
for  Cf/Cftat  each  streamwise  location  for  each  solution. 

With  the  two-point,  windward  scheme  leading  edge  error  influ¬ 
ences  the  solution  far  downstream.  The  three-point,  windward 
scheme  damps  the  error  out  more  effectively,  converging  on 
better  yalues.  The  three-point,  windward  scheme  is  therefore 
more  accurate.  Part  of  the  reason  for  reduced  leading  edge 
error  in  the  three-point,  windward  scheme  is  this  scheme  uses 
three  analytic  points  at  the  leading  edge.  The  two-point, 
differencing  method  only  uses  two  analytic  points  at  the 
leading  edge.  In  the  analytic  solution,  the  Yg  metric  is 
known  exactly.  It  is  calculated  by  knowing  the  exact'  form 
of  the  bound a r y- 1  ay e r  shape.  Eq  (48)  computes  the  Yg  instead 
of  calculating  it  from  numerical  differencing.  This  takes 
the  numerical  error  out  of  the  Y^  metric.  With  a  completely 
analytic  Yg  metric  and  density  and  viscosity  forced  to  remain 
constant,  the  best  Cf  ratio  obtained  with  a  61x30  grid  is 
.99361.  Figs.  10  compares  this  completely  analytical  case  to 
the  case  with  Y«  computed  numerically  except  at  the  first 
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three,  streamwi.se  locations.  Density  and  viscosity  are  held 
constant  in  both  runs.  The  completely  analytic  case  is 
definitely  more  accurate  for  predicting  flow  in  the  boundary 
layer.  However,  if  the  analytic  metric  is  used,  the  boun¬ 
dary-layer  code  could  only  be  used  with  a  solution  grid  which 
has  been  scaled  according  to  a  square-root  of  X  function  as 
in  Eq  (38).  This  scaled  grid  is  accurate  for  similar  bounda¬ 
ry  layers.  Since  the  boundary-layer  code  needs  to  be  applied 
to  more  general  problems,  the  numerical  metric  is  used.  Also 
plotted  in  Fig.  10  is  the  solution  with  the  numerical  7^  and 
changing  density  and  viscosity.  There  is  a  very  alight  error 
between  the  two  numerical  Yg  solutions  caused  by  allowing 
density  and  viscosity  to  vary.  This  error  does  not  affect 
the  results,  significantly.  The  new  method  of  determining 
and  three-point  differencing  on  Y^  significantly  improves  the 
solutions  computed  by  the  bounds ry- 1  aye r  code.  The  grid 
dependence  of  the  finite-difference  solution  also  cannot 
be  ignored. 


Grid  Dependence . 

The  number  of  points  in  the  solution  grid  and  the  way 
they  are  generated , affect  the  boundary- 1 ayer  solution.  As 
the  number  of  points  in  the  solution  grid  increases,  the 
solution's  acOuracy  increases,  and  the  computer  time,  or 
iterations,  required  for  a  converged  solution  also  increas¬ 
es.  Fig.  11  shows  the  solutions  for  three,  different,  size 
grids.  By  increasing  the  number  of  grid  points  Cf/Cft  ap~ 
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proaches  one.  However,  the  61x39  grid  requires  75  tine 
steps  and  206. 78  CPU  seconds  to  converge  on  the  solution. 

The  61x30  grid  takes  20  tiae  steps  and  29  CPU  seconds.  The 
21x19  grid  takes  5.97  CPU  seconds  for  20  time  steps.  In¬ 
creasing  the  number  of  grid  points  does  increase  accuracy 
but  with  a  considerable  increase  in  computer  time. ' 

Using  the  correct  guess  for  a  solution  grid  also  in¬ 
creases  the  accuracy  of  a  non-adaptive  grid  solution.  Lange 
uses  an  optimized  grid  developed  by  Hodge,  and  the  non-adap¬ 
tive-grid  solutions  in  this  study  use  an  optimized,  power-law 
grid  to  solve  the  boundary-layer  code  (10).  However,  the 
adaptive  grid  solutions  in  this  study  use  an  exponential  grid 
generated  by  a  solution  of  the  grid  equation  in  the  normal 
direction.  Pig.  12  demonstrates  the  grid  dependence  of  the 
boundary- 1  aye r  code.  The  exponential  grids  are  generated 
with  Q  equal  to  -.3.  The  61x30  power  law  grid  results  are 
much  better  than  the  results  for  either  a  61x30  or  a  21x19 
exponential  grid.  The  type  of  grid  chosen  for  the  boundary- 
layer  solution  does  affect  the  accuracy  of  the  solution.  The 
power  law  grid  is  accurate  for  b ound a ry- 1  aye r  problems  since 
the  boundary  layer  is  close  to  a  power-law  shape  except  at 
the  surface.  But  the  exponential  grid  has  the  control  func¬ 
tion,  Q,  which  can  be  varied  depending  on  the  characteristics 
of  the  flow  solution.  Optimizing  Q  yields  an  exponential  grid 
which  gives  accurate  results  and  does  not  require  any  fore- 
knowlege  of  the  exact  solution. 
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Non- Ad i  ab  a  t  i  c .  Wall  Effects. 

The  test  cases  previously  mentioned  use  input  conditions 
which  result  in  a  nearly,  adiabatic  flow.  The  wall  tempera¬ 
ture  is  close  to  the  adiabatic,  wall  temperature  calculated 
from  Eq  (28).  The  approximately,  adiabatic  wall  magnifies 
the  effect  of  roundoff  error  for  the  calculation  of  h.  The 
finite-difference  solution  uses  Eq  (55)  to  find  h.  The 
denominator  in  Eq  (55)  is  very  small  for  a  nearly,  adiabatic 

wall.  Dividing  q  by  a  very,  small  number  produces  good  re- 

/ 

suits  which  give  an  h/href  near  1.0,  but  computer,  round-off 
error  causes  the  computed  h  to  overshoot  the  values.  Fig.  13 
compares  h/hrej  for  the  nearly,  adiabatic  wall  case  and  a 
non-adiabat ic  case.  The  non-adiabat ic  case  uses  a  wall 
temperature  of  550°R,  instead  of  530°R  as  in  the  adiabatic- 
wall  case.  The  two  solutions  use  a  61x59  solution  grid  and 
an  analytic  definition  of  T^.  The  non- ad i ab a t  ic  case  does 
not  overshoot  the  theoretical  solution,  but  it  does  have 
some  error  due  to  compressibility. 

The  non- ad i ab a t i c  wall  case  highlights  the  behavior  of 
the  finite-difference  solution  toward  h.  With  the  non-adia- 
batic  variation  in  wall  temperature,  compressibility  effects 
are  much  greater  than  previous  cases.  The  density  ranges 
from  .9963075  to  .99941725  across  the  boundary  layer.  There¬ 
fore,  there  is  some  error  in  the  solution  due  to  this  com¬ 
pressibility.  Fig.  14  shows  Cf/Cf^  is  not  affected  much  by 
the  wall  temperature  change.  However,  h  is  affected.  The 
non- ad  i  ab a t  ic  cases  also  male  the  solution  for  the  total 
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Pig.  14  Non-adlabattc  Wall  Effecta,  Cf/Cf^. 


enthalpy  much  more  difficult.  For  the  adiabatic  cases,  it 
takes  20  iterations  to  converge  to  EPS  of  lxlO-^,  and  U  was 
the  usually  the  property  that  determined  DELE.  For  the 
non-adiabat ic  cases,  it  takes  over  80  time  iterations  to 
converge  to  EPS  of  1x10”^.  The  finite-difference  solution 
definitely  has  more  trouble  resolving  H  for  the  non-adiabat ic 
cases.  The  finite-difference  solution  does  converge  to 
values  close  to  the  theoretical  values,  though.  Using  the 
same  finite-difference  solution  with  s  Thomas- a  1 gor i thm, 
iteration  method.  Hodge  has  computed  ratios  for  Cf/Cft  and 
h/hf%f  of  .9998  and  1.00,  respectively  (8).  Therefore,  the 
finite-difference  solution  of  the  incompressible,  boundary- 
layer  solution  closely  reproduces  Blasius',  exact  solution. 

Optimized  Boundary  Laver  Code . 

The  adaptive-grid  solution  of  the  s t e ady- s t a t e ,  boun¬ 
dary-layer  problem  finds  the  solution  grid  which  best  re¬ 
solves  boundary-layer  flow.  The  finite-difference  solution 
of  the  b ounda ry- 1  ay e r  equations  is  very,  grid  dependent.  If 
an  optimized,  solution  grid  is  found,  the  accuracy  of  the 
computed  solution  increases.  Using  Powell's  method  to  mini¬ 
mize  the  sum  of  the  squares  of  tie  strearawise  velocity  gradi¬ 
ents  optimizes  the  grid  and  increases  the  accuracy  of  the 
computed  solution.  The  increase  in  accuracy  is  more  easily 
seen  with  a  coarse  grid,  so  the  computer  solutions  use  a 
21x19  exponential  grid.  After  showing  that  Powell's  optimi¬ 
zation  method  produces  an  optimized  grid,  this  study  test* 
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the  nee  of  different  minimizing  functions  in  the  optimize- 

7.  7.  9 

tion.  These  functions  ere  the  suia  of  U  ,  0^  ,  or  • 

Minimizing  some  functions  produces  better  optimizetion  then 
others.  The  number  of  iteretions.  KT,  end  other  input  pern- 
meters  such  es  the  initiel  guess  for  Q  also  effect  the  opti¬ 
mizetion.  If  the  initiel  guess  for  Q  is  too  fer  off.  Pow¬ 
ell's  method  may  not  reech  e  converged  solution  in  e  smell 
number  of  iteretions.  Without  e  converged  solution,  the 
optimized  solution  for  the  psrticular  function  may  not  be 
reached.  Powell's  method  does  minimize  most  of  the  input 
functions  in  very  few,  iteration  steps. 

Optimizetion  Performance . 

The  computer  solution  generates  .  the  solution  grid  using 
tii  grid  equation  in  Eq  (35).  The  optimizetion  finds  the 
control  function,  Q,  which  generates  a  grid  that  minimizes 
the  desired  function  and  increases  the  accuracy  of  the  boun¬ 
dary-layer  solution.  The  flexibility  of  Powell's  method 

2  2  2 

allows  the  sum  of  U„  ,  0„_  ,  or  to  be  minimized. 

q  qq  qqq 

Comparison  of  these  three  cases  determines  which  minimization 
gives  the  best  results.  As  in  the  non-adaptive  case,  a 
better  solution  comes  closer  to  a  value  of  1.0  for  Cf/Cft. 

In  addition,  the  optimization  which  reduces  the  error  in  the 
computed  solution  is  a  better  method.  The  root-mean-square 
(RMS)  of  the  difference  between  the  computed  solution  and  the 
von  Karman-Pohlhausen  approximate  solution  measures  the  error 
in  the  computed  solution.  The  different  optimization  solu¬ 
tions  are  also  compared  to  the  Blasius',  exact  solution. 
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With  these  comparisons,  it  is  possible  to  discover  how  the 
adaptive  grid  is  best  used  to  get  more,  accnrate,  competed 
solutions  to  general,  boundary-layer-type  problems. 

Table  III  summarizes  test  cases  of  the  adaptive-grid 
solution  tor  the  incompressible  flow  over  a  flat  plate.  The 
input  parameters  for  solution  of  the  bounds ry- 1  aye r  code 
within  the  optimization  do  not  change.  The  only  change  in 
the  bounda ry- 1  ay e r  code  inputs  is  the  solution  grid  calculat¬ 
ed  in  the  optimization.  The  boundary- layer  code,  input 
parameters  are  20  time  steps,  4  iterations  per  time  step,  a 
time  step  of  1x10®  seconds,  error  tolerance  of  1x10“®,  a  Pr 
of  .72,  and  recalculation  of  the  SOR  parameter  every  fifth 
time  step.  These  are  the  same  inputs  as  run  5  of  the  non- a- 
daptive,  incompressible  solution.  Optimizing  21x19  and 
61x30  exponential  grids  is  tested.  Setting  the  minimum  and 
maximum  values  of  Y  at  the  x  location  where  the  grid  equation 
is  solved  determines  the  physical  size  of  grid.  These  Y 
values  are  YMIN  and  YMAX  in  Table  III A.  The  grid  is  then 
scaled  in  the  streamwise  direction  with  the  Y  spacing  solved 
by  the  grid  equation  at  the  x  location  chosen.  For  this 
study,  the  grid  equation  is  solved  at  the  trailing  edge 
of  the  surface.  The  parameter,  Q,  shows  the  initial  input 
for  the  control  function.  The  results  of  using  a  non— adap¬ 
tive,  exponential  grid  with  the  input  Q's  are  given  in  Table 
I V .  KT  represents  the  number  of  iterations  allowed  to  find 
each  X.r  value.  PCT  is  the  percentage  of  Xr  added  and  sub¬ 
tracted  to  Xr  to  get  X+  and  X“ .  The  error  tolerance  between 
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TABLE  I I IB 

Adaptive  Exponential  Grid  Solutions-  Results 


Run 

Cf/Cft 

(leading  edge) 

Cf/Cft 

(trailing  edge) 

F 

( initial-final) 

RMS(U-Uexact) 

( initial-final) 

1 

.99135 

.98781 

.093644-. 028557 

.029654-. 028736 

2 

.99135 

.98781 

.093644-. P2255 

.029654-. 02836 

3 

.99135 

.98781 

.093644-. 028557 

.029654-. 028736 

4 

.99589 

.99590 

.220564-. U00611 

.032002-. 033511 

5 

.98473 

.98394 

72.308  -  70.934 

.019075-. 019598 

6 

.98624 

.98281 

2.2295  -  2.2073 

.029654-. 029671 

7 

.98636 

.98292 

?.2i;44-2.::.07 

.029654-. 029671 

8 

.99039 

.98637 

2.0612  -2.0403 

.029797-. 029946 

9 

.99107 

.98760 

.3C723  -.23333 

.029654-. 029764 

10 

.99053 

.98707 

1.16154-1  0384 

.029654-. 029925 

11 

.99139 

.98735 

.24103  -.10459 

.029654-. 028933 

12 

.99084 

.98' 38 

.66444  -.56026 

.029654-. 029861 

13 

.99138 

.98784 

.16734  -.063586 

.029654-. 028809 

14 

.64983 

.64920 

72.308  -71.503 

.019075-. 019381 

15 

.65324 

.65259 

72.308  -70.934 

.019075-. 019598 

16 

.94844 

.94325 

6.176  -  6.1275 

.02993  -.029903 
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TABLE  I I IB  (cont.) 


Cf/Cf. 


RMS(O-Dexact) 

initial- 


027888-^02377 


UU 


.99135 

.96920 

.97324 

.95975 


09365- . 022259 

.02 9654-. 028732 

24102- . 10455 

.029654-. 028834 

.229  -  2.203 

,029654-. 026815 

093654-. 022544 

.0,29654-.  02  .>7  36 

.98781 

.96591 

.96990 

.95656 


093655- . 022544 
866760- . 605740 
690502- . 690086 
,38803-3.73679 


.029654-.028736 
.036949-. 037135 
.036949-. 036962 
.036949- .038336 


These  ere  set  et  lxlO-5.  The  lest  input  parameters  shown  in 
Teble  III  ere  the  weighting  peremeters,  WH1 ,  WH2  ,  end  W 03 . 
These  weightings  determine  which  velocity  derivative,  or 
combinstion  of  velocity  derivetives,  is  the  minimized  func¬ 
tion.  The  minimized  function  is  represented  by  P  end  is 
defined  by 

F  «  2+WH2*£D  2+fH3*£U  *  (66) 

r\  HH 

TH1  is  the  weighting  on  minimizing  JCU^2.  WH2  is  the  weight¬ 
ing  on  minimizing  ZU^2,  end  VB3  is  the  weighting  on  minimiz¬ 
ing  £lJ„  2  The  velues  in  Teble  IIIB  ere  the  renge  of  re- 
Hi 

suits  showing  the  performence  of  the  particular,  edeptive- 
grid  solution. 

Powell's  method  minimizes  the  input  function.  As  e 

result,  the  computed,  boundery  leyer  solution  gets  closer  to 

theoreticel  velues.  This  study  seeks  to  minimize  the  coapu- 

tetionel  error  by  minimizing  truncetion  error.  From  en 

enelysis  of  the  truncetion  error  terms,  minimizing  0„__ 

B'l'l 

should  reduce  solution  error.  Therefore,  this  study  concen- 
tretes  on  minimizing  the  input  function,  »  with  Pow¬ 

ell's  minimization  method.  Solution  runs  1-4,  shown  in 
Table  III,  summarize  the  results  of  optimizing  the  grid  with 
The  differences  between  the  runs  ere  the  various 
input  parameters.  All  of  the  solutions  show  e  decrease  in 
.  represented  by  F,  from  the  start  of  the  program  to 
its  converged  solution.  Each  of  the  solutions  also  comes 
closer  to  the  theoretical  solution  than  the  initial  grid. 


The  range  of  values  for  C  f  /  C  f  t  is  closer  to  1.0  than  similar 
results  with  either  the  initial,  exponential  grid  or  the 
non-adapt  ive ,  power-law  grid  solution  shown  in  Table  IV. 
Additionally,  the  optimized  solatiou  reduces  the  overall 
error  in  0.  In  run  1,  the  RMS  of  U-Uex#ct  across  the  whole 
domain  is  .029654  at  the  start  of  the  minimization.  The 
final  RMS  ®”®exaot  value  is  .028736.  This  is  a  3.096%  de¬ 
crease  in  the  0  error.  Fig.  15  shows  how  close  the  solution 
comes  to  the  exact  u  velocity  profile  over  a  flat  plate. 
Therefore,  the  Powell's  minimization  method  applied  to  the 
input  function  reduces  the  error  in  U  and  increases 

nnn 

the  accuracy  of  the  computed  solution  by  minimizing  the  input 
function,  F . 

Although  minimizing  meets  the  objective  of  re¬ 

ducing  the  error  of  the  computed  solution  and  increasing 
its  accuracy,  it  may  not  be  the  only  function  which  can  meet 
these  objectives.’  Using  ZU_^,  ZU__^,  or  some  combination  of 

n  ’in 

the  three  velocity  gradients  as  minimizing  functions  produces 
different  results  than  those  optimizations  on  EU___^,  Each 

nnn 

of  the  different  cases  performs  similar  to  minimizing  EU^.^ 
by.  decreasing  F  and  increasing  the  solution  accuracy.  Howev¬ 
er,  the  error  in  U  is  not  decreased  in  all  cases.  Minimizing 
does  net  decrease  the  error  in  U.  Optimizing  this 
function  increases  the  error  in  U  by  .06%.  Cf/Cft  for  this 
solution,  run  6,  is  also  below  that  of  the  runs  1-4.  Since 
the  function  did  decrease,  an  increase  in  the  error  is  not 
expected.  Closer  examination  of  the  solution  shows  the 


optimization  is  not  totally  converged.  The  solution  was 
therefore  rerun  with  a  higher  number  of  iterations.  A  Q 
closer  to  the  results  of  previous  runs  was  also  input  to  help 
the  solution  converge  faster.  This  solution,  run  8  in  Table 
III,  did  give  better  resnlts  for  Cf/Cft.  The  results  are 
still  not  as  good  as  those  in  run  1.  The  error  in  D  also 
increased  by  .496%.  Minimizing  ZO^2  does  keep  more  points  in 
the  boundary  layer  than  other  solutions.  Sun  1  keeps  11 
points  in  the  boundary  layer.  Suns  6-8  keep  14  points  in  the 
boundary  layer.  With  more  points  in  the  boundary  layer,  the 
solution  should  resolve  flow  conditions  better  and  produce 
better  results.  However,  this  does  not  occur.  Fig.  16  shows 
the  velocity  profile  produced  by  minimizing  ZO^2.  It  is 
close  to  the  theoretical  answer.  However,  Fig.  17  shows,  the 
velocity  profile  solved  by  minimizing  ZD^2  is  farther  from 
the  exact  solution  than  either  minimizing  ZB„„  2  or  ZU_  2 . 
Regardless  of  the  number  of  points  in  the  boundary  layer,  the 
solution  from  optimizing  ZO^2  is  not  better  than  optimizing 
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Optimizing  with  ZU„  2  is  better  than  optimizing  5"lu _ 2 

in  some  ways.  Minimizing  ZU^2  in  run  11  gets  almost  as 
close  to  the  theoretical  results  as  runs  1-4.  Cf/Cft  starts 
out  higher  than  case  1  but  recovers  to  values  lower  than 
run  1.  The  decrease  in  error  for  D  is  only  2.766%.  Minimi¬ 
zing  resolves  the  leading  edge  gradients  better  ,o  get 

better  values  at  the  leading  edge.  But  it  does  not  produce 
better  results  at  the  trailing  edge  where  the  leading  edge 


76 


•t 


errcr  does  not  affect  the  solution  as  much.  Over  the  entire 
dontin,  the  sum  of  U_  2  does  not  reduce  the  error  in  D  as 
well  as  minimizing  the  sum  of  U^,^2.  Pig.  17  shows  that 
minimizing  minimizes  the  error  in  0  better  and  pro¬ 

duces  a  veloc i ty  prof  ile  closer  to  the  exact  solution  than 

minimizing  either  the  sum  of  U_  *  or  the  sum  of  U„  . 

•  n’t  ’l 

Using  each  of  the  three,  different,  velocity  gradients 
separately  as  the  minimizing  function  in  Powell’s  method 
shows  the  sum  of  U_„_  to  be  the  best  minimizing  function. 
Using  combinations  of  the  three  velocity  gradients  also  does 
not  improve  the  solution.  In  runs  6-13,  different  weightings 
are  applied  to  the  different  gradients.  The  weighting  deter¬ 
mines  how  much  a  part  of  the  minimized  function  the  specific 
gradient  is.  When  ^U^2  is  any  part  of  P,  the  velocity  error 
increases  and  the  computed  solutions  are  worse  than  run  1 
where  only  )U„_  2  is  minimized.  Combining  the  sum  of  U  2 
and  U  ,  as  in  run  13,  gives  slightly  better  results  for 
C  f  /  Cf  j. .  In  the  combination,  the  ^U^  2  contribution  resolves 
the  leading  edge  error  while  ^U^.^2  raises  the  trailing  edge 
values.  The  combination  of  the  two  functions  thus  produces  a 
grid  which  gets  the  closest  values  to  the  Blasius',  theoreti¬ 
cal  solution.  The  error  in  U  only  decreases  by  2.85%. 
Minimizing  >U__  2  still  reduces  the  error  in  U  better  than 
any  of  the  functions  tested.  There  may  be  some  combination 
of  the  ]>Unl12  and  Jun„n2  not  tested  which  may  reduce  the  error 
in  U  better  in  addition  to  computing  better  results  for  Cj. 

Several  of  the  input  parameters  affect  the  behavior  of 
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the  optimization  method.  Inputs  KT,  PCT,  and  Q  affect  the 
convergence  of  the  method,  and  YMAX  affects  the  accuracy  of 
the  solution.  Increasing  KT  only  affects  those  optimizations 
which  need  many  iterations  to  converge.  Powell's  method  hat 
many  levels  at  which  it  iterates.  The  iteration  of  and 
iteration  of  the  total  solution  are  two  of  the  levels.  Al¬ 
though  a  problem  may  fail  to  converge  at  one  level,  succeed¬ 
ing  levels  take  out  this  error  to  converge  with  a  small  KT. 
rising  a  large  KT  for  quickly  converging  problems  converges 
the  solution  in  the  iteration  of  Xr.  All  steps  which  follow 
merely  check  this  convergence.  This  checking  uses  computer 
time  but  adds  nothing  to  the  solution.  The  latter  steps  of 
the  method  are  therefore  wasted.  This  occurs  in  run  2. 
However,  for  slow  converging  problems,  increasing  KT  is 
essential  to  reach  Converged  answers.  Optimizing  on 

M 

requires  many  iterations.  For  example,  run  6  does  not  con¬ 
verge  with  a  low  KT.  Runs  7  and  8  converge  with  only  double 
the  KT  of  run  6.  The  results  of  the  converged  problem  are 
better  than  those  of  the  unconverged  problem.  In  run  8,  the 
initial  Q  is  changed  to  help  convergence.  If  the  initial  Q 
is  very  far  away  from  the  final,  optimum  Q,  the  solution 
needs  more  iterations  to  converge.  Runs  14  and  IS  are  solu¬ 
tions  with  an  initial  Q  of  .3.  From  solutions  1-3,  the 
optimum  Q  for  the  same  minimized  function  and  grid  size  is 
-.144.  In  runs  14  and  15,  Q  only  decreases  from  .3  to  .297 
in  three  iterations.  This  solution  is  not  converged.  The 
solution  grid  resulting  is  not  optimized,  as  shown  by  the 
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values  of  Cf/Cf^.  When  the  initial  Q  is  changed  to  .05  in 
run  17,  the  solution  converges  to  -.144  the  same  as  runs 
1-3.  This  shows  positive  guesses  for  Q  can  converge,  but  the 
Initial  guess  must  be  relatively  close  to  the  optimum  Q  for  a 
small  KT.  Large,  negative  guesses  for  Q  also  do  not  con¬ 
verge.  W4*h  Q  initially  set  at  -.6,  the  solution  converges 
on  Q  of  -.594.  This  is  a  long  way  from  the  run  1  solution  of 
-.144.  Negative,  initial  guesses  for  Qcannot  be  too  far. 
away  from  the  final  Q.  There  is  a  another  limit  on  negative, 
Q  values.  Q  cannot  at  any  time  in  the  iterations  be  allowed 
to  get  too  small  for  proper  resolution  of  D  in  the  finite 
difference.  In  the  21x19  grid,  this  *3  is  -.8.  For  a  61x59 
grid,  Q  of  -.3  compresses  the  grid  too  much.  The  U  veloci¬ 
ties  for  a  very,  compressed  grid  are  virtually  the  same  for 
the  y  locations  next  to  the  wall.  Therefore,  the  velocity 
gradients  go  to  zero,  and  the  solution  diverges.  Large, 
negative  guesses  for  Q  should  be  chosen  carefully.  Choosing 
a  large  PCT  also  can  cause  the  solution  to  diverge.  PCT 
determines  how  far  away  from  X,  X  +  and  X-  are.  This  optimi¬ 
zation  method  performs  like  any  iterative  scheme.  If  guessed 
values  are  too  widely  spaced,  more  iterations  are  required  to 
converge  to  a  solution  ,  or  the  solution  may  even  diverge. 

Run  3  which  uses  PCT  of  .5  converges  on  the  same  answer  as 
run  1 .  Run  1  converges  in  6  iterations,  but  run  3  converges 
in  9.  The  larger  PCT  is  therefore  not  as  efficient,  unless 
the  initial  guess  is  very  far  away  from  the  optimum  Q. 

Choosing  PCT,  KT,  and  the  initial  Q  correctly  leads  to  a 
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converged  answer,  but  the  size  of  the  physical  grid  affects 
the  accuracy  of  the  solution.  Choosing  YMIN  and  YMAX  deter¬ 
mines  the  physical  size  of  the  grid.  The  exponential  grid  is 
built  within  these  limits.  For  most  of  the  solutions,  YMAX 
is  .015.  This  is  approximately  46  for  Re  of  2x10^.  In  run 
4,  YMAX  is  dtvreastd  to  .004,  or  approximately  5.  ,  This  puts 
the  majority  of  the  points  inside  the  boundary  layer.  The 
optimization  solves  the  boundary  layer  better  when  only  the 
bou.dairy  layer  is  included  in  the  domain.  As  expected, 
solution  4  is  very  close  to  the  Blasius'  solution.  If  YMAX 
is  expanded  to  8<J  as  in  run  5  ,  the  solution  is  not  as  good 
as  solutions  with  thinner,  solution  grids.  Therefore,  the 
choice  of  YMAX  does  affect  the  optimized  solution.  If  solu¬ 
tions  are  made  w  i  t  h  v  ar  i  sib  1  e  ,  flow  conditions,  YMAX  should  be 
varied  to  keep  the  same  48  grid  size  in  the  y  direction. 
Otherwise,  the  results  are  not  compatible. 

Various,  input  parameters  and  minimizing  different  func¬ 
tions  affects  the  optimized  solution.  Changing  the  flow 
conditions  or  the  number  of  grid  points  does  not  change  the 
behavior  of  the  solution  for  giver  minimizing  functions.  In 
runs  13-22,  putting  in  differe:  Reynold's  numbers  for  cases 

minimizing  ,  or  tests  the  effect  of  differ¬ 

ent,  flow  conditions.  The  solutions  from  runs  18,  19,  20,  or 

21  are  no  different  at  P.e  of  5x10^  than  the  corresponding 

c  v  2 

runs  1,  6,  11,  or  2  Re  of  2x10°.  For  minimizing  ' 

there  is  no  difference  for  any  of  the  solutions  in  runs  1,  2, 

I*.  21,  or  22.  Changing  Reynold's  number  has  no  affect  on 
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which  minimizing  function  optimizes  the  solution  best. 

Changing  the  number  of  grid  points  alters  the  accuracy 
of  the  solution.  In  run  26,  a  61x30  grid  definitely  gives 
better  values  of  Cf/Cft  than  the  21x19  grid  used  in  run  1. 
This  is  expected.  Decreasing  the  number  of  grid  points 
decreases  the  accuracy  of  the  solution.  Runs  23-23  are  two 
to  three  percent  less  accurate  than  similar  runs  1,  7,  and  11 
with  the  21x19  grid.  However,  using  the  courser,  21x10  grid 
changes  the  relative  accuracy  of  the  three  minimizations 
procedures.  Minimizing  the  sum  of  U  *  gives  better  results 
for  the  course  grid  than  does  minimizing  the  sum  of  or 

The  optimum  value  of  Q  also  changes  with  the  number 
of  grid  points.  A  larger  number  of  grid  points  requires  an 
optimum  Q  closer  to  0.  The  21x10  grid  in  run  23  converged  on 
a  Q  of  -.33361.  The  21x19  grid  in  run  1  converged  on  a  Q  of 
-.14426,  while  the  61x30  grid  in  run  26  required  an  optimum  Q 
of  -.0848.  As  the  number  of  grid  points  increases,  the 
optimum  solution  comes  closer  to  a  uniform  grid.  Therefore, 
the  number  of  grid  points  does  affect  the  accuracy  and  the 
method  which  best  optimizes  the  b ound a ry- 1  ay e r  solution. 

Powell's  method  optimizes  the  incompressible,  boundary- 
layer  solution.  The  sum  of  U  A  is  the  best  minimizing 
function  to  minimize  the  error  in  the  streamwise  velocity. 
Parameters,  YMAX,  KT,  PCT,  and  the  initial  Q  determine  how 
fast  and  what  values  the  solution  converges  on.  Changing  the 
Reynold's  number  of  the  flow  problem  does  not  affect  the 
behavior  of  the  optimization.  Increasing  the  number  of  grid 
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points  does  increase  the  accuracy  of  the  solution  and  does 
change  the  optimization  behavior  for  course,  solution  grids. 
Incompressible,  flow  problems  can  be  solved  more  accurately 
with  an  adaptive  grid  using  fewer,  grid  points.  The  optimi¬ 
zation  can  also  be  applied  to  resolve  more  variable,  compres¬ 
sible,  or  turbulent,  bound  a r y- 1  ay e r  problems. 


t  B 


Comp  r  e  s  s  ib  1  e  Flow 

The  behavior  of  th*  optimization  and  boundary- 1  aye r 
codes  has  already  been  examined  for  the  incompressible,  flow 
problems.  Since  the  programs  were  initially  set  up  to  handle 
compressible,  flow  problems  and  had  to  be  fooled  to  treat 
incompressible  flow,  the  extension  of  the  optimization  to 
supersonic  and  hypersonic,  compressible,  flow  problems  is  an 
easy  step.  The  two,  largest  problems  are  defining  an  initial 
grid  for  the  f in i t e- d i f f e r en c e  solution  and  comparing  the 
results  with  previous,  theoretical  data.  In  the  incompressi¬ 
ble  cases,  the  boundary- 1 ayer  thickness  determined  the  physi¬ 
cal  size  of  the  computational  domnin.  TMAX,  the  height  of 
the  grid  at  the  surface’s  trailing  edge,  was  set  at  approxi¬ 
mately  46  to  get  consistent  results.  For  the  compressible 
cases,  the  boundary  layer  is  not  constant.  The  physical  size 
of  the  boundary  layer  changes  with  flow  velocity.  Also,  the 
transformed,  bound  a  r  y- 1  a'-e  r  thickness  used  with  the  compres¬ 
sibility  transformation  changes  with  the  ratio  Tw/Te,  not 
with  the  velocity.  Therefore,  as  the  compressible  solutions 
are  applied  to  different  velocities  and  temperature  condi- 
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tions,  some  method  of  defining  a  consistent,  transformed- 
grid  size  mnst  be  found.  The  othei  problem  with  the  compres¬ 
sible  cases  is  determining  the  accuracy  of  the  computed 
solution.  There  is  no  exact  solution  for  compressible  flow 
over  a  flat  plate  or  wedge.  There  is  no  equivalent  to  the 
incompressible,  Blasins’  solutions.  Eckert  theory  predicts 
approximate  values  for  Cf  and  h  across  a  high-speed  flow. 
Also,  Van  Driest’s  study  of  compressible,  bounds ry- 1  ay e r  flow 
over  a  flat  plate  provides  a  theoretical  development  to 
compare  with  the  computed  results.  These  comparisons  show 
the  validity  and  limitations  of  the  boundary-layer  code. 

The  major  difference  between  the  compressible  and  incom¬ 
pressible,  boundary- layer  solutions  is  the  activation  of 
the  compressibility  t r ans f o rma  i  i on .  As  flow  velocity  in¬ 
creases  and  temperature  variations  occur  across  the  boundary 
layer,  compressibility  affects  the  flow  more  and  cannot  be 
ignored.  This  study  assumes  compressibility  effects  to  be 
important  beginning  at  about  Mach  .05.  Above  Mach  .05,  the 
Dorodnitsyn  transformation  applies  to  the  problem.  The 
problem  is  set  up  with  the  same  equations  as  the  incompressi¬ 
ble  problem.  However,  a  new,  b ound a ry-  1  ay e r  thickness,  A(X), 
is  calculated  with  Eq  (16).  Above  the  transformed,  boundary 
layer,  the  flow  is  inviscid.  Initial  conditions  are  the  same 
as  the  incompressible  case,  except  that  the  initial  thermal 
b ound a r y- 1  ay e r  thickness  is  assumed  equal  to  the  transformed 
bo und a r y- 1  ay e r  thickness.  Since  A(X)  varies  depending  on 
the  edge  conditions,  differing  Mach  number,  Reynold's  num- 


ber,  or  temperature  conditions  require  a  variable  definition 
on  the  physical  size  of  the  solution  grid.  In  the  incompres¬ 
sible  cases  studied,  YMAX,  the  height  of  the  solution  grid  at 
the  surface's  trailing  edge,  is  a  fixed  input.  For  the 
compressible  cases,  YMAX  varies  with  the  size  of  A(X).  It 
can  be  set  at  any  multiple  of  A(X)  depending  on  the  require¬ 
ments  of  the  flow  problem.  The  effect  of  using  different 
multiples  of  A(X)  on  the  solution  is  shown  in  Fig..  18.  This 
figure  shows  the  computed,  velocity  profile  at  the  trailing 
edge  for  different,  solution  methods.  Similar  to  the  incom¬ 
pressible  cases,  the  closer  YMAX  is  to  A (X)  the  better  the 
solution  is  to  theoretical  results.  The  solution  for  YMAX  at 
2A(X)  is  closer  to  the  theoretical  result  than  3A(X)  or  4A,(X) 
solutions.  The  application  of  the  Dorodnitsyn  transformation 
leads  to  computed  solution  which  behaves  like  the  incompres¬ 
sible  solutions  in  other  ways,  also. 

Although  the  compressible,  boundary-layer  problem  does 
not  have  an  exact  solution  similar  to  the  Blasius ’  solution 
for.  flat  plate  flow,  this  study  compares  the  computed  solu¬ 
tion  to  theoretical  results  obtained  by  E.  R.  Van  Driest. 

Van  Driest  solves  the  laminar,  compressible,  boundary- 1  ay e r 
equations  for  flow  problem  over  a  flat  plate  using  Crocco's 
method  (24:1).  Van  Driest' !  conditions  on  the  flow  problem 
are 

Prandtl  number  was  taken  at  0.75,  the  specific 
heat  constant,  and  the  Sutherland's  law  of  viscos¬ 
ity-temperature  variation  was  assumed  to  represent 
viscosity  data  starting  with  an  ambient  temperature 
of  -67.6°F  (24:  1)  . 
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Pig.  18  Expanded  Mach  4  Velocity  Prof lle-ru^  minimized 


Van  Driest':  assumptions  coincide  with  approximations  made 
in  this  study.  Van  Driest  solves  the  flow  problem  for  many, 
different,  flow  parameters.  Rather  than  show  a  comparison 
with  the  extensive  results  of  Van  Driest's  development,  this 
stndy  examines  a  few  cases  to  compare  computed  solutions  to 
Van  Driest's,  theoretical  solutions. 

The  computed  solutions  are  very  close  to  Van  Driest's 
solutions.  The  cases  solved  are  for  Pr  of  .75,  the  ratio  of 
Sutherland  constant,  198. $°R,  to  free  stream  temperature  at 
.  505  ,  and  Tw/T00  at  1.  This  study  tests  supersonic  and  hyper¬ 
sonic  cases  for  Mach  2  to  8.  Pigs.  19  and  20  compare  the 
computed  velocity  and  temperature  profiles  across  the  plate 
for  the  Mach  range.  These  solutions  set  the  upper  limit  of 
the  domain  at  3A(X)  away  from  the  plate  and  optimize  on  the 
sum  of  •  The  computed,  velocity  profiles  are  very  close 

to  Van  Driest's  profiles.  In  fact,  the  Mach  4  theoretical 
velocity  profile  lies  exactly  over  the  computed  solution. 

Fig.  21  shows  several,  different,  optimized  solutions  for 
Mach  4.  All  of  the  optimizations  are  so  close  that  an  expan¬ 
sion  of  the  graph  is  needed  to  see  any  difference  in  the 
solutions.  Fig.  18  shows  the  Mach  4  velocity  profiles  that 
result  from  optimizing  on  the  sum  of  and  on  the  sum  of 

Tltin 

2  2 
U  .  Optimizing  with  the  sum  of  is  slightly  better 

2 

than  optimizing  with  the  sum  of  U.._  .  This  is  consistent 

with  the  results  of  the  incompressible  cases.  Although  the 
Mach  8,  velocity  profile  shows  more  error  than  the  Mach  >4 
case,  the  Mach  8,  velocity  profile  is  also  close.  Since  this 
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study  is  restricted  to  flows  where  a  perfect-gas  assumption 
is  valid,  some  variation  from  the  Van  Driest's,  theoretical 
model  is  expected  for  Mach  numbers  above  6.  The  tempera¬ 
ture  profiles  in  Fig.  20  are  not  quite  as  exact.  The  com¬ 
puted.  temperature  solutions  are  less  than  the  theoretical 
solution  at  some  points.  Since  the  boundury-layer  code 
reduces  the  error  in  0  by  minimizing  the  sum  of  ,  the 

computed  solution  should  be  close  to  the  theoretical  solu¬ 
tion.  However,  the  temperature  profile  results  show  the  code 
does  not  necessarily  produce  optimized  temperature  results 
for  Prandtl  number  not  equal  to  one.  It  does  achieve  the 
goal  of  reducing  the  error  in  IT. 

The  computed  results  also  come  close  to  values  of  h 
predicted  by  high-speed  Eckert  theory.  For  high-speed  prob¬ 
lems,  Eckert  recommends  using  an  approximate  temperature, 

T*.  to  calculate  an  approximate  density,  ,  and  viscosity, 
ji  ,  for  the  compressible,  boundary  layer.  Eq  (33)  gives  a 
value  for  T* .  The  approximate,  heat  convection  coefficient, 
h*,  is  then  calculated  by  putting  the  *  conditions  in  Eq  (32) 
a  «  an  approximation  for  h  across  the  compressible,  boun¬ 
dary  layer.  Eq  (32)  only  applies  to  lminar  flow  with  Re  not 
greater  than  5x10^  (12:213).  By  restricting  flow  conditions 
to  Re  less  than  this  limit,  the  computed  results  for  h  come 
close  to  Eckert's,  high-speed  predictions.  Fig.  22  shows  the 
computed  solutions  for  Mach  2,  4,  6,  and  8  relative  to 
Eckert's  solution.  A  solution  at  Mach  14.24  over  a  20°  wedge 
is  also  included  to  show  the  solution  does  work  over  a  wedge. 
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Both  the  theoretical  and  confuted  solution  are  divided  by 

lijefi  the  heat  convection  coefficient  calculated  for  free- 

stream  ccnditicns.  All  of  the  solutions  shown  in  Pig.  22 

* 

optimize  the  solution  grid  with  the  sum  of  U  over  3^(1). 

If  a  1A(X)  solution  grid  is  used,  the  computed  h  should  get 
closer  to  Eckert's  theory.  Fig.  23  shows  this  is  not  the 
case.  Fig.'  23  presents  h/hrej  results  for  Mach  4  with  YMAX 
at  different  multiples  of  A(X).  The  increase  in  h/hre£  for 
1A(X)  solutions  results  because  less  than  the  true,  trans¬ 
form  d,  bounds ry- 1  ay e r  thickness  is  included  in  the  domain. 
The  .nitiul  guess  for  the  boundary- layer  thickness,  A(X),  is 
too  small.  The  computed  solution  actually  uses  about  1.3A(X) 
as  .he  point  where  inviscid  flow  occurs  and  D  e / U  ^  n  £  is  one. 
Com: arisons  of  the  results  in  Fig.  23  indicate  that  the 
closer  YMAX  is  to  the  true,  boundary  layer,  the  better  the 
results  for  h  are.  However,  the  entire  boundary  layer  must 
be  included.  Since  the  initial  guess  is  not  precise,  it  is 
diff.crlt  to  pick  the  YMAX  position  to  get  only  the  boundary 
lav  r.  Including  too  much  of  the  inviscid  flow  above  the 
bound !  ry  layer  also  dilutes  the  boundary- 1 ayer  results. 

Some  problems  may  require  the  solution  )f  more  than  the  flow 
inside  the  boundary  layer.  When  the  problem  is  set  up,  the 
desired  accuracy  and  size  of  the  problem  domain  must  be  con¬ 
sidered  for  the  best  results  possible.  The  optimized,  boun¬ 
dary-layer  solution  provides  the  most  flexibility  for  dealing 
with  compressible,  laminar,  boundary- layer,  flow  problems  and 
still  achieving  results  close  to  theoretical  predictions. 


Chanter  VI :  Cone  lus  ions 

In  incompressible,  flow  problems,  the  theoretical  solu¬ 
tion  gives  a  cine  that  the  solution  grid  should  be  a  power- 
law  grid.  For  most  flow  problems,  however,  theoretical  solu¬ 
tions  do  not  exist  to  give  a  guess  for  the  solution  grid. 
Being  able  to  input  a  general,  exponential  grid,  which  can  be 
optimized  to  give  the  best  flow  solution,  allows  the  most 
general  application  of  a  b ound a r y- 1  ay e r  code. 

The  bounds ry- 1  ay e r  code  presented  in  this  study  is  an 
improvement  over  Lange's,  b ound a r y- 1  ay e r  code.  For  the  non- 
adaptive  solutions,  the  modifications  made  to  Lange's  code 
improve  the  solutions.  Non-adaptive  solutions  for  Cj  calcu¬ 
lated  by  the  new  code  are  closer  to  the  Blasius’,  exact  solu¬ 
tion  for  flow  over  a  flat  plate.  In  addition,  the  new  code 
includes  a  minimization  method  which  optimizes  an  adaptive 
grid  to  find  better  solutions  of  boundary-layer  problems. 

Optimizing  the  adaptive  grid  with  Powell's  method  pro¬ 
duces  more  accurate,  boundary-layer  solutions,  but  the  inputs 
must  be  chosen  carefully  to  get  the  best  optimization.  The 
bounds ry- 1  ay e r  solution  is  very,  grid  dependent.  The  number 
of  grid  points,  physical  size,  and  type  of  grid  definitely 
affects  the  results  of  the  finite-difference  solution.  In 
both  non-adaptive  and  adaptive  cases,  the  solution  with  the 
largest  number  of  grid  points  gives  the  best  results.  The 
problem  whose  domain  includes  less  of  the  inviscid  region 
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outside  the  boundary  layer  produces  the  better  results.  For 
example,  problems  that  included  only  16  in  the  normal  direc¬ 
tion  had  better  results  than  problems  that  used  46  as  the 
height  of  the  grid.  More  grid  points  are  in  the  boundary 
layer  with  only  16  as  the  domain  limit.  In  the  non-adaptive 
cases,  the  p owe r- 1  aw- t yp e ,  solution  grid  gives  better  results 
than  an  arbitrary,  exponential  grid.  However,  the  power-law 
grid  only  applies  to  a  specific,  bound a r y- 1  ay e r ,  problem 
type.  Although  the  initial,  exponential  grid  input  may  give 
worse  results  than  the  power-law  grid,  using  Powell's  method 
to  find  the  optimized,  exponential  grid  produces  the  best 
results.  The  initial,  control  function,  Q,  must  be  chosen 
close  to  the  final,  optimized  value.  If  it  is  chosen  too  far 
away  from  the  final  Q,  Powell's  method  will  not  converge. 
Convergence  is  also  affected  by  the  number  of  iterations 
allowed  and  the  spread  between  minimization  guesses.  If  too, 
few  iterations  are  done,  the  solution  does  not  converge  on  an 
optimized  value.  If  too  many  iterations  are  used,  all  of  the 
convergence  takes  place  in  the  initial  steps  of  Powell's 
method.  This  wastes  all  subsequent  calculations  which  check 
the  convergence  and  is  inefficient.  Choosing  a  large  spread 
between  minimization  guesses  also  is  less  efficient.  Pow¬ 
ell's  method  calculates  a  central,  functional  value,  f,  two 
functional  values  on  either  side  of  the  central  value,  f+  and 
f  ,  and  the  minimum  functional  value  of  a  quadratic  fit 
through  the  previous  three  functional  values,  fB.  If  the 
spread  between  the  f,  f+,  and  f~  is  too  large,  the  solution 
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takes  store  iterations  to  converge.  The  successive  iterations 
tend  to  continually  overshoot  the  final  value  until  getting 
very  close.  A  large  spread  nay  be  useful  in  problems  where  a 
good  initial  guess  is  not  known,  though.  The  function  mini¬ 
mized  also  affects  the  efficiency  of  the  adaptive-grid  solu¬ 
tion. 

Powell's  method  can  minimize  any  input  function.  This 

y 

study  optimizes  the  solution  grid  with  the  sum  of  D 

,  and  of  combinations  of  these  three  functions.  Of 

2 

the  first  three  .functions,  minimizing  the  sum  of  0^^^  best 
minimizes  the  error  in  D .  It  produces  better  results  for 
Cf.  The  solution  reduces  the  RMS  of  the  error  between  the 
computed  U  and  the  von  Ka r man-Poh 1 hau s en  approximation  of  D 
most.  The  computed  velocity  profile  is  also  closest  to  Bias- 

y 

ins',  exact,  velocity  profile.  Minimizing  the,  sum  of 
does  not  produce  good  results.  The  minimization  does  not 
converge  on  an  optimum  grid  with  this  function.  It  does  put 
more  points  into  the  boundary  layer,  but  the  results  are 
worse  for  Cf  and  the  velocity  profile.  In  fact,  minimiz- 
ing  ^  ,  increases  the  error  between  the  computed  U  and  the 
von  Ka rma n-Po h 1 h a u s e n  solution.  If  ^  is  included  with  any 
of  the  other  derivative  functions,  the  results  are  always 
worse.  This  is  not  the  case  with  minimizing  Mini¬ 
mizing  with  ^  gets  results  very  close  to  the  results  of 

minimizing  the  sum  of  Minimizing  better  re- 

solves  the  flow  at  the  leading  edge,  while  minimizing 
better  resolves  flow  at  the  trailing  edge.  Consequently,  the 
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combination  of  5n  and  shows  improved  results  over 

u  nn  l  nnn 

minimizing  separately.  The  ^  still  produces  the 

least  error  in  U,  however.  Some  other  combination  of  the 
Kn2  *nd  tte  may  produce  better  results.  Changing 

Reynold's  number  for  the  problem  does  not  significantly 
change  which  minimized  function  produces  the  best  results.' 
Changing  the  number  of  grid  points  in  the  y  direction  seems 
to  make  the  ^  the  best,  minimizing  function.  For  prob¬ 

lems  with  smaller  numbers  of  grid  points,  the  affect  of  mini- 
mizing  the  increases.  The  number  of  grid  points  also 

affects  the  optimized  Q  value  and  the  initial  Q  guess.  Grids 
with  more  grid  points  converge  on  Q's  closer  to  0.  Since  the 
optimization  inputs  do  affect  the  accuracy  and  behavior  of 
the  minimization  process  and  final  solution,  the  problem  must 
be  examined  carefully  to  determine  the  inputs  which  give  the 
best,  adaptive-grid  solution. 

Using  the  ad  ap  t  i  ve- g.r  i  d  solutions  in  compressible,  boun¬ 
dary-layer  problems  also  produces  good  results.  Computed 
velocity  profiles  are  an  excellent  match  to  Van  Driest's, 
theoretical  solutions  for  high-speed  flow  over  flat  plates. 
The  temperature,  profile  results  are  not  quite  as  exact,  but 
the  computed,  temperature  profiles  follow  Van  Driest's,  theo¬ 
retical  profiles  closely.  The  a d ap t iv e- gr id  results  also 
compare  well  with  high-speed  Eckert  theory  predictions  for 
the  heat  convection  coefficient.  Like  the  incompressible 
cases,  minimizing  the  produces  the  best  results  in 

compressible  cases.  Throughout  the  range  from  Mach  2  to 
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Mach  8,  the  compnted  h  does  not  vary  from  Eckert's  theoreti¬ 
cal  valnes  by  more  than  three  percent.  This  is  good  correla¬ 
tion  to  an  approximation  like  Eckert  theory.  Therefore,  the 
adaptive  grid  solutions  sncceed  in  minimizing  the  error  in  U 
to  give  excellent  velocity  profiles  for  high  speed,  compres¬ 
sible  flow  and  good  agreement  with  high  speed  Eckert  theory. 

The  bounds ry- 1  ay e r  code  presented  in  this  study  is  an 
improvement  over  previous  methods.  The  non-adapt ive .  grid 
solutions  are  better  than  Lange’s  solutions  and  are  close  to 
Blasius.'  exact, incoapressible  solution.  Adaptive-grid  solu¬ 
tions  improve  on  the  solutions  using  a  non-adaptive  grid  and 
use  fewer  grid  points.  The  adaptive  grid  also  provides  the 
f 1  ex ib i 1 ity  1  to  optimize  the  grid  for  any  desired  function  and 
for  general  boundary- 1 ayer  problems  that  may  not  have  exact 
solution  such  as  high-speed  compressible  or  turbulent  flow. 
The  adaptive  grid  successfully  solves  compressible,  flow 
problems.  The  application  of  the  optimized,  exponential, 
adaptive-grid  method  to  these  compressible  problems  shows 
that  the  method  can  be  applied  to  more  complicated  problems 
than  incompressible  flow  over  a  flat  plate.  Using  a  non- 
adaptive  grid  structured  around  a  grid  optimized  for  the  flow 
solution,  limits  the  application  of  the  boundary- 1  aye r  meth¬ 
od.  The  adaptive  grid  has  no  such  limitations.  Further 
application  of  the  ad ap t i v e- g r  i  d  method  could  solve  even  more 
complicated  problems.  The  solution  of  these  problems  may 
extend  the  engineer's  knowledge  of  hypersonic,  flow  problems. 
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Chanter  VII:  Recommendat ions 


This  study  has  not  investigated  all  possible  applica¬ 
tions  or  affects  of  the  adaptive-grid  optimization.  Further 
studies  of  the  optimization  procedure  could  make  some  changes 
which  could  improve  the  efficiency  of  the  optimization. 

Other  changes  to  the  bound a r y- 1  ay e r  code  could  also  expand 
the  applicability  of  the  code  to  a  much  larger  class  of  prob¬ 
lems. 

Two  improvements  to  the  optimization  code  could  improve 
its  efficiency.  The  biggest  shortcoming  in  the  adaptive-grid 
solution  is  that  it  requires  much  more  computer  time  than  a 
non-adaptive  scheme  since  a  solution  of  the  boundary- 1  aye r 
code  is  required  for  each  iteration  of  the  method.  Even 
though  the  number  of  grid  points  may  be  less,  the  number  of 
iterations  required  for  the  adaptive-grid  solution  increases 
the  overall,  computer  time.  If  a  Thoma  s a  1  go  r  i  t  hm,  iterative 
solution  is  used  instead  of  the  SOR,  the  time  required  to 
compute  the  bound a r y- 1  ay e r  solution  decreases.  Chen  has 
included  a  quad-diagonal  solver  in  this  same,  boundary-layer 
code  to  study  turbulent,  flow  problems  (5).  This  would  dra¬ 
matically  reduce  the  computer  time  needed  for  adaptive-grid 
solutions.  Also,  minimizing  some  other  combination  of  the 
sum  of  0^  and  the  sum  of  U  may  produce  better  overall 

results  than  cases  tested  in  this  study.  Other  minimizing 
functions  could  also  be  tested.  The  parameters  that  form  the 
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control  function  could  also  be  changed.  This  study  splits  Q 
into  two  parts  and  minimizes  each  part.  Q  could  be  split 
into  more  parts,  or  a  linear  functional  relationship  between 
the  parts  could  be  used  to  form  Q.  Many  different  parame¬ 
ters  can  be  varied  to  find  the  best  optimization  method.  The 
flexibility  of  Powell's  method  provides  at  any  avenues  which 
might  result  in  better  answers  for  various  problems. 

There  are  additions  that  can  be  made  to  apply  the  opti¬ 
mization  grid  to  a  larger  class  of  problems.  To  further 
reduce  the  error  in  the  boundary  layer  code  and  make  it  a 
true,  two-dimensional  problem,  the  adaptive  grid  could  be 
applied  to  both  streamwise  and  normal  directions.  This  study 
only  optimized  in  the  normal  direction  and  then  scaled  the 
exponential  stretching  in  the  streamwise  direction.  Turbu¬ 
lence  models,  and  non i $ o the rma 1  wall  effects  could  be  in¬ 
serted.  The  ax i symme t r ic ,  bound a ry- 1  ay e r  equations  could  be 
used  to  solve  flows  over  axisymmetric  bodies.  This  would 
allow  the  solution  of  a  wide  variety  of  bounda ry- 1  ay e r  flows 
with  the  adaptive  grid.  The  possibilities  for  using  the 
basic  optimization  method  presented  in  this  study  are  limit¬ 
less.  But  much  more  study  is  required  to  determine  the  ap¬ 
plication  of  ad ap t  ive- gr  id  and  boundary-layer  solutions  to 
investigations  o f  hypersonic,  flow  regimes. 
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The  bound* ry- 1  aye r  equations  under  most  conditions  are 
density  dependent.  Comp r e s s ib 1  i ty  makes  the  solution  of 
the  non-linear,  boundary- 1  aye r  equations  more  difficult.  If 
the  equations  are  trausformed  into  equations  similar  to  the 
incompressible,  boundary- layer  equations,  a  solution  might  be 
easier.  To  transform  the  compressible  equations,  some  form 
of  a  compressibility  transformation  is  used.  The  transforma¬ 
tion  puts  the  density  dependence  into  the  varibles  of  the 
equations.  In  the  physcial  space,  the  dimensional  variables 
are  x',  y',  and  t*.  The  variables  are  non-d  imens  ional  ized 
using 


The  transformed  variables  are 

ry 

X  =  x,  T  =  J  o  *  dy,  t '  =  t  (A- 2) 

0  P'e 

The  non-dimensional,  compressible,  axisymmetric  boundary- 
layer  equations  are 

Continuity:  rm9_£.  +  9  r^o  u  +  9  r  mp  v  =  0.  (A-3) 

9 1  9  x  9  y 

9  ou  +  9  ou^  +  9  puv  =  - 9_p  +  9 / _ji  9_a 

9t  9x  9y  9x  9y\Re  9 y 


Momentum : 


(A-4a) 


Momentum : 


3n  =  0 

i'y 


( A-4b ) 


Energy:  3  o  H  +  3  ouH  +  3  o  vH  *>  3_p  +  3  /  n  3  H 

3t  3x  3y  3t  dy^RePr  3y 

+  _3 

®y 

where  b  =0  for  two-dimen* ional  problems  and  m=l  for  axi- 
symmetric  problems  (Appendix  B) .  The  transformed  variables 
mast  be  snbstitated  into  the  compressible,  boundary-layer 
equations.  To  substitute  for  the  partial  derivative  terms 
use  the  following  relationships: 

35=3W91  +  3J<LX*  3J  *  3J  31  +  O  £Y,  3W  =  3W  (A-6) 

3 x  3 X  3 x  3 Y  3 x  3y  3X  3y  3Y  3y  3t  3t 

where  3_X=  3_Y=  3 1  =  3t=  3X=0  and  3_Y=  p 
3t  3t3x3y3y  3y 

W  in  the  above  equations  represents  any  variable  of  interest. 
After  expanding  Eq  ( A- 3 )  using  the  chain  rule  and  then  sub-1 
st  i  tut ing  relations  from  Eq  (A-6),  the  transformed,  continu¬ 
ity  equation  becomes 

rB3_£+  rm3j>  3Y+ttr“  3£+prB/3D+iS  3l\  +  r®v/<LfiL  £1=0 

3 1  3 Y  3 1  3 x  \3X  3Y  3x/  \3Y  3y/  3Y  3y 

(A-7) 

The  desired  form  of  the  incompressible,  ax i symme t r i c ,  conti¬ 
nuity  equation  is 

3r®P  +  3  rBV  =  0  (A-8) 

3X  3  Y 

The  continuity  equation  in  ax i s ymme t r ic ,  physical  space  co¬ 
ordinates  is 


RePr 


3  u  /  2 ) 

3y  J 


( A-5 ) 
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V"  v 


a  -n'i  +  a  rmv  =  o 
a  x  a  y 

Compare  Eqs  (A-8)  and  (A-9)  to  determine  what  0  is, 

$ rm_  =  r"3u  3X  =  rm3u  =  rBaP 

ax  ax  dx  ax  ax 

Therefore, 


(A-9) 


( A- 1 0  ) 


( A-l 1 ) 


Taking  oat  the  dD/3X  term  in  Eq  (A-7),  the  remaining  terms 
are  the  second  term  in  Eq  (A-8) . 


rB3j£+urB  a_£+prmajJ  3Y  +  3  rmo  v  =  9  rttv  =  p  r 111  B V 
dt  dx  3Y  dx  dy  dy  dT 


where  p=  3Y/3y 


( A-12 ) 


rMa  *  Y  +  r*uaL_Y  +  rBaa  aYi8r"cv  -  prBaV  (A-13) 
3y3t  ByBx  dy  8x  By  By 


rmp v=  rBV  -  rBaY/at  -  urBaY/3x 


{ A-14 ) 


There  fore , 


V=  pv+  BY/Bt  +  uBY/Bx 


( A- 1 5 ) 


The  transformed  velocity  components,  U  and  V,  satisfy  Eq 
(A-8).  When  they  are  substituted  into  Eqs  (A- 4)  and  (A- 5), 
equations  similar  to  the  incompressible,  momentum  and  energy 


equations  result  (15:81-83).  These  equations  are 


Moaentut 


au+cao+vau  =  - 1_  aj>+a (o u  a  p\ 

at  ax  aY  p  ax  by  ^Re  by) 


a  p  =  o 
aY 


Energy:  3 R+P3  H+V9H 

at  ax  a y 


( A- 1 6  a ) 


(A- 16b ) 


0=1  3p +a  (  DU  aj\ 

3 Y  p  at  aY\RePr  BY  ) 

+  9  /pu(Pr-l)  3P!l2j\  (A-17) 

9Y  \  RePr  3Y  / 


107 


The  two-dimensional  analysis  of  boundary- 1  aye r  equations 
is  easily  applied  to  axisymmetric  shapes.  Axisymmetric 
shapes  have  one  axis  such  that  a  plane  passed  perpendicular 
through  this  axis  slices  out  a  circle  of  varying  radius. 

Fig.  B-l  shows  a  general,  axisymmetic.  coordinate  system. 

Ary  location  on  the  surface  is  described  by  coordinates  x.  y, 
and  ♦.  The  scale  factors  are  h^.  Velocity  components  are 
Vi .  hj  and  Vj  are 

hj  =  1  +  y/Rj  Vx  -  u 

h2  =  1  V2  -  v  (B-l) 

hj  =  R0(x)  +  ycosa  V3  =  0 

Rq  is  the  local  radius  of  the  surface  and  varies  with  x.  Rj 

is  the  local  radius  of  curvature.  For  a  cone,  a  is  0°,  and 
there  is  no  x-dependence  in  y.  The  only  non-unity,  scale 
factor  is  I13  which  becomes  R  (x),  or  just  r.  To  make  the 
two-dimensional  equations  more  applicable  to  two-dimensional 
and  conical,  axisymmetric  problems,  the  superscipt  m  is  add¬ 
ed.  For  two-dimensional  problems,  m  is  0.  For  axisymmetric 
problems,  m  is  1.  In  this  study,  the  two-dimensional,  incom¬ 
pressible,  boundary-layer  equations  were  transformed  into  a 
computational  plane  where 

$  =  ?(X.  Y ,  t) 

q  =  q  (  X  ,  Y ,  t )  ,  ( B-2 ) 

T  =  t 
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After  the  metric  relationships  are  applied,  the  transformed, 
two-dimensional,  incompressible,  boundary-layer  equations  are 

Continuity:  3  iaP  5r  +  3  rmU  rw  +  9  rmV  tty  =  ®  (B-3) 

a?  35  35 

Momentum:  Ut  +  U(  0$$x  +  Cq,lx)  +VDq,'Y  “  1  /  Re  t pp  ( ^  R  Y  <  B—  4  ) 
Energy:  Ht  +  +  +  VHq,'Y  *  ( B-  5 ) 

pt/p  +  (  ou  H  nY)  ’'Y  +  t^U-LEilJLL  ny  ^  ^  n  ^  RY 

Pr Re  n  2RePr  2 

The  momentum  and  energy  equations  are  not  affected  by  the 
changing  radius  of  a  cone's  surface.  Consequently,  the  fi¬ 
nite-difference  equations  for  finding  U  and  H  do  not  change. 
The  solution  of  V  requires  integrating  the  continuity  equa¬ 
tion  which  does  depend  on  the  changing  radins. 

The  continuity  equation  is  rearranged  to  isolate  V. 

The  ,  r e sul t  ing  eqaution  is  then  integrated  to  solve  for  V. 

The  integrable  equation  is 

<rmV)n  =  ^_(rn0)5  +  Y£_(rmU)n  (B-6) 

x5  x? 

The  right-hand,  side  term  with  (rmII)^  is  differenced  with  a 
three-point,  windward  scheme  about  the  j+1/2  point.  The 
second  term  splits  up  into  parts. 

Y.  (rmU)  =  1_[  Y,  ( rmH>  -  f(  rmU)(Y«)  dr,  ]  J  (B-7) 

X’  X?  J  j-1 

The  integral  is  evaluated  using  a  trapezoidal  rule  with  Y 
averaged  about  the  j-1/2  point  prior  to  the  integration. 
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Y5/X$(rmD)n  =  l/X^  {t(Y5rmD)  i#  j.j]  ,  (B-8) 

-  .5t(r«U)iij  +  (r'Il)iij.1]lY5iij  "YU.j-l]} 

-  1  /  ( 2Xg  )  f(Y^rmD)  lf  j-(Y5rBU)  i#  (B-9) 

+  Y$i,  j-l(rmn)  i.  j“Y$i.  j(r"D)  i.  j"l] 

-  1/lx««<Yti.j^Y5ifJ-1)t(.-n)i#r(r-n)i(<5iJj) 


With  the  other  terms  of  added  in,  the  finite  difference 

for  V.  .  is 
1 »  J 


vi.j  -  i/«*">i.4  {<-"v>i,j-i 


+  l/(ax4)C.S(Yi#J-IlfJ.1)C-3Kr»1»1<J  +  (r-0)liJ.1l 
+4I(r«D)i.lij  +  (r»D)..1(j.iH(,»D)i.2>j  +  (r‘U)i.2|j.1]J 
+  <Y?i.j+T$i.j-l)nr“U)i.j-(raiJ)i.j-l1)}  (B"11> 


This  finite  difference  is  then  iterated  to  solve  for  V  at 
each  location  in  the  domain. 
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Appendix  D:  Optimization  Program 
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PROGRAM  MNTRER 

C  THIS  PROGRAM  USES  POWELL'S  METHOD  OF  CONJUGATE  DIRECTIONS 
C  TO  MIN  THE  TRUNCATION  ERROR  IN  A  BOUNDARY  LAYER 
C  THIS  PROGRAM  WILL  COMPUTE  THE  FOLLOWING: 

C  NUMERICAL  GRID  SOLUTION  USING  POWELL'S  METHOD  AND 
C  THE  CONTROL  FUNCTION  P  TO  MINIMIZE  TRUNCATION  ERROR 
C  IN  THE  TRANSFORMED  ETA  DIRECTION  OF  A  BOUNDARY  LAYER 
C  GRID 

C  IN  ALL  CASES  THE  THIRD  DERIVATIVE  (DUUU)  IS  CALCULATED  NUMERICALLY 

C 

C 

DIMENSION  U£XACT(21 .0:20) ,P(19) ,X(21) ,DUUU(2l) 

DIMENSION  Y ( 21 . 19) , YPLUSDL (21 , 19) , YMINDL (21 , 19) 

DIMENSION  P3( 19) ,D3UC21) ,Y3(21,19) ,U3<21,0:20) 

DIMENSION  AAA (19) ,BBB(19) , CCC(19> , DDD(19) , GGG<19) , WWW (19) 
DIMENSION  PPLUSDL(19) ,PMINDL(19) 

DIMENSION  U ( 21 , 0 : 20 ) , UPLUSDL (21,0:20), UMINDL (21,0:20) 

DIMENSION  A (5) , ASAVE(S) ,A3(5) ,APLUSDL(5) , AMINDL(S) 

DIMENSION  ABAR(S) ,21(2,2) , DELTA (21) ,ABAR2(5) 

DIMENSION  POLD( 19) 

0PEN(UNIT=8,FILE*'RESTART' ) 

OPEN ( UNIT=9, FILE* 'WALLQ' ) 

OPEN ( UNIT=10 , FILE* ' FIELD' ) 

CPEN(UNIT*11,FILE*'GRID') 

0PEN(UNIT=12,FILE='HREF' ) 

OPEN ( UNIT* 13, FILE* 'CFCF' ) 

REWIND  6 
REWIND  8 
REWIND  9 
REWIND  10 
REWIND  11 
REWIND  12 
REWIND  13 
C 

CALL  SECOND (CPI) 

READ  (7, • )  IMAX , JMAX , KMAX , ERMXL 
READ  (7,*)  YMIN,YMAX, PERCENT 
READ  (7, •)  NMAX.ITMAX 
READ  (7, •)  ERPOWMX , ICOL , I RST 
READ  (7,*)  (A(IR) ,IR=1,NMAX) 

READ  (7,*)  WH0,WH2,WH3 
READ  (7,«)  2 
CALL  DATE(ADATE) 

CALL  TIME(ATIME) 

WRITE(6, 100)  ADATE, ATIME 
100  FORMAT ( 'l',/,5X,A10,2X,A10) 

WRITE(6,200)  IMAX, JMAX, YMIN.YMAX, KMAX, ERMXL, PERCENT, 

1  (A(IR),IR=1,NMAX) 

200  F0RKAT< IX, 'NUMBER  OF  GRID  PTS,  IMAX*' ,I3,3X, ' JMAX=' , 

1  13,/, IX, 'MIN  Y  VALUE,  YNIN*' ,E13.5,/,1X, 

2  'MAX  Y  VALUE,  YMAX*' ,E13.5,/,1X, 

3  'MAX  NUMBER  OF  ITERATIONS  FOR  ONE', 

4  'PARAMETER  OPTIMIZATION,  KMAX*' ,15,/, IX, 

5  'MAX  PERCENT  ERROR  IN  ONE  PARAMETER', 

6  '  OPTIMIZATION  FOR  CONVERGENCE,  ERMXL*', 

7  E13.5,/, IX, 'PERCENT  CHANGE  IN  PARAMETERS  USED', 

8  'TO  CALCULATE  QUAD  EQ  FOR  ONE  PARAMETER  OPT*', 

9  E13. 5,/, IX, 'INITIAL  A  ' ,5E13.5,/,/) 

WRITE  (6,209)  WH0,WH2,WH3 

209  FORMAT  (IX, 'THE  WEIGHTING  FOR  DU  *  ',E7.2,2X, 

1  'DUU  =' ,E7.2,2X, 'DUUU  *',E7.2,/) 

WRITE  (6,211)  ITMAX 

211  FORMATdX, 'MAX  NUMBER  OF  ITERATIONS, IN  OUTER  LOOP', 
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ooon  oo 


1  ■  'TO  POWELL  METHOD,  ITMAX=',I5) 

213  FORMATUX,  'OPTIMIZED  COLUMN  IS',1X,I3> 

PI=AC0S(-1.) 

JMMWMAX-1. 

SETTING  UP  BOUNDARY  CONDITIONS 
XSTEP«1/(FL0AT(IMAX)-1.) 

DO  5  1=1, IMAX 
X ( I > “FLOAT  < I - 1 ) »XSTEP 
CONTINUE 
Y(IC0L,1)*YMIM 
Y  < ICOL , JMAX ) *  YMAX 

THIS  SECTION  CALCULATES  THE  OPTIMIZED  NUMERICAL  SOLUTION 
USING  POWELL'S  METHOD  (1964) 


DO  50  J*1 , NMAX 
ASAVE( J) *A( J) 

DO  70  IR=1,NMAX 
ZI(J,IR)*0.0 
CONTINUE 
2I(J, J)3A(J)»0.1 
IF  (ZI(J.J) .EQ.O.O)  ZI ( J, J) 3 .01 
CONTINUE 
IB=0 

CALL  PCAL ( A , NMAX , P , JMAX  > 

CALL  EXPGRD ( X , Y , P , A AA , CCC , DDD , GGG , WWW , YMIN , 

L  YMAX,IMAX, JMAX, ICOL) 

CALL  BLIMP ( X , Y , IMAX , JMAX , U , IB , UE , DELTA , RE ) 

WRITE (6,*)  'DELTA3  ', DELTA (IMAX) 

IF  (YMAX.LT.Z*DELTA(IMAX) )  THEN 
YMAX=DELTA(IMAX)*Z 
Y ( ICOL , JMAX ) = YMAX 
GO  TO  55 
ENDIF 

CALL  NWP  (U,DUUU, Ml, IMAX, JMAX, ICOL, WHO, WH2.WH3) 

H=H1 

CALL  POHL ( Y . DELTA , UE , I MAX , JM AX , UEX ACT ) 

WRITE  <6  402) 

FORMAT  (3X,'I',4X,'J',09X,'P',12X.'X',12X.'Y',10X, 

L  'UEXACT' , 8X, 'U' ,11X, 'U-UEXACT' ) 

DO  611  I?1 , IMAX 
IF<I. EG. 2. OR. I. EG. IMAX)  THEN 
DO  612  J=1 , JMAX 

WRITE  (6,502)  I,J,P(J),X(I),Y(I,J) , UEXACT ( I , J) , 

1  U(I,J),U(I,J) -UEXACT ( I , J ) 

FORMAT  (1X.2I3.6E13.5) 

CONTINUE 
ENDIF 
CONTINUE 
WRITE  (6,602) 

FORMAT  (IX,'  P  IS  JUST  A  FUNCTION  OF  THE  INITIAL  A', 

1  'VALUES  AND  NOT  OF  RLAMBA  OR  ZI',/) 

CALL  ERROR(U, UEXACT, IMAX, JMAX) 

ITRT=1 

IF  (ITRT.GT.NMAX)  GO  TO  580 
DO  370  ITERATE=ITRT, ITMAX 
HOLD3H 
IRSAVEsO 
HDELMAX=0.0 
ERR0RPW=0.0 
DO  150  IR=1 ,NMAX 

ERRORPW=ERRORPW+ABS <  A ( IR ) - ASAVE ( IR) ) 
ASAVE(IR)=A(IR> 

CONTINUE 

ERRORPW=ERRORPW/NMAX 

IF  (ERRORPW.LE.ERPOWMX. AND. ITERATE. NE.l)  GO  TO  580 
WRITE  (6,122)  ITERATE 
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122  FORMAT  <1 X,/,/,/,' . BELOW  ITERATE®'  ,15, 

1  '  . ') 

IF(IRST.EQ.l)  THEN 
C 

C  STEP  1  OF  POWELL'S  ALGORITHM:  FOR  IR*1 ,2, . . .NMAX  CALCULATE  RLAMBA 
C  SO  THAT  H(A+RLAMBA»2I)  IS  A  MINIMUM  AND  DEFINE  A(IR) =A( IR-1) *RLAMBA»2I 
C 

IR1-1 

DO  ISO  IR=IR1,NMAX 
RLAMBA=0.0 
DLTA=PERCENT»RLAMBA 
IF  (ABS(DLTA) .LE. .05)  DLTA*.05 
RLAMPL»RLAMBA*DLTA 
RLAMMN=RLAMBA-DLTA 
C 

C  ’  , 

C  ITERATION  FOR  ONE  PARAMETER  OPTIMIZATION 
C 

DO  250  K*1 ,KMAX 
RLAMOLD* RLAMBA 
DO  ISO  J  =  1  ,NMAX 

APLUSDH J) »A< J) ♦RLAMPL*ZI ( J , IR) 

AMINDL(J) *A(  J)-»RLAMMN«ZI ( J, IR) 

ABAR( J) SA( J) +RLAMBA»2I ( J , IR ) 

180  CONTINUE 

CALL  PCAL( APLUSDL, NMAX, PPLUSDL, JMAX) 

CALL  PCAL(AMINDL, NMAX, PMINDL, JMAX) 

CALL  PCAL ( ABAR, NMAX, P, JMAX) 

CALL  EXPGRD ( X ,  Y ,  P , AAA , CCC , DDD , GGG , WWW , YMIN , 

1  YMAX  IMAX  JMAX  ICOL) 

CALL  EXPGRD  (X, YPLUSDL.PPLUSDL, AAA, CCC, DDD, GGG, WWW, YHIN, YMAX, 
1  v  IMAX, JMAX, ICOL) 

CALL  EXPGRD  (X. YMINDL.PMINDL, AAA, CCC, DDD, GGG, WWW, YMIN, YMAX, 

1  IMAX, JMAX, ICOL) 

C 

CALL  BLIMP < X , Y , IMAX , JMAX , U , ITERATE , UE , DELTA , RE ) 

CALL  BLIMP(X,YPLUSDL, IMAX, JMAX, UPLUSDL, ITERATE, UE, DELTA, RE) 

CALL  BLIMP < X , YMINDL , IMAX , JM AX , UMINDL , ITERATE , UE , DELTA , RE ) 

CALL  NEWP1  (U,P, UPLUSDL, UMINDL, DUUU,RLAMPL,RLAMMN, 

1  RLAMBA, RMSDUUU,H,K,PERCENT,IMAX, JMAX, 

2  ABAR2, NMAX, X.Y, AAA, CCC.DDD, GGG, WWW, YMIN, YMAX, ICOL, 

3  ITERATE .  UE , DELTA , P3 , Y3 , D3U , U3 , A , 21 , IR , WHO , WH2 , WH3 ) 

WRITE(6,800)  IR,(A(J),J=1,NMAX) 

800  FORMAT (IX, 'ABOVE  RESULTS  ARE  FOR  INNER  OPT  AND', 

1  'A(IR)  PARAMETER,  IR=',I5,/,'  THE  OLD  A(I>', 

2  'VALUES  ARE=',5E13.5> 

WRITE  (6,801)  (2I(J,IR),J=1,NMAX) 

801  FORMAT  (IX.'  THE  OLD  ZI(J,IR)  VALUES®' ,5E13. 5,/, 

1  IX,'  POLD  IS  CALCULATED  BY  USING', 

2  '  A® A (OLD) ♦ (RLAMBA  OLD)»(ZI  OLD)',/,/) 

ERRORL® ABS( RLAMOLD-RLAMBA ) • 100 . 0 

IF  (RLAMOLD . NE . 0 . 0 )  ERRORL® ERRORL/ ABS(RLAMOLD) 

IF  ( ERRORL. LT.ERMXL)  THEN 
WRITE(6,»)  'K-CONVERGED' 

GO  TO  290 
ENDIF 

250  CONTINUE 

290  CONTINUE 

DO  310  J*1,NMAX 
A( J) =A( J) ♦RLAMBA«ZI ( J , IR) 

310  CONTINUE 

C 

C  STEP  2:  FIND  THE  INTEGER  IRSAVE  SO  THAT  H(A(IR-1))- 
C  H( A(IR) )  IS  A  MAXIMUM,  AND  DEFINE  THE  MAXIMUM  AS  HDELTNAX 
C 

HDELTA*HOLD-H 

HOLD=H 

.  IF(HDELTA.GT.HDELMAX)  THEN 
HDELMAX®HDELTA 
IRSAVE® IR 
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END  IF 
CONTINUE 


C 

C  STEP  3*.  CALCULATE  H3=H(2*A(NMAX)-A(0) ) ,  AND  DEFINE  H1=H(A<0)). 

C  AND  H2CH(A(NMAX) ) 

C 

163  H2*H 

IRST=0 

DO  320  IR=1,NMAX 
A3(IR) *2«A(IR) -ASAVE(IR) 

320  CONTINUE 

CALL  PCAL<A3,NMAX,P3,JMAX) 

CALL  EXPGRD(X,Y3,P3, AAA, CCC, DDD, GGG. WWW, YMIN, YMAX, IMAX, JMAX, ICOL) 
CALL  BLIMP(X,Y3, IMAX, JMAX , U3, ITERATE, UE, DELTA, RE) 

WRITE < 6, • )  'THIS  NWP  CALLED  IN  STEP  3' 

CALL  NUP ( U3 ,  D3U ,  H3 ,  INAX , JMAX ,  ICOL , WHO ,  W*i2 ,  WH3 )  C 
C  STEP  4:  IF  EITHER  H3.GE.H1  AND/OR 

C  (2»<H1-2«H2*H3) • (H1-H2-HDELMAX) *«2) .GE. (HDELMAX* (H1-H3) *»2) 

C  USE  THE  OLD  DIRECTIONS  21  FOR  THE  NEXT  ITERATION  AND 

C  USE  A(NMAX)  FOR  THE  NEXT  A<0)  (GO  TO  STEP  1),  OTHERWISE 

C  GO  TO  STEP  5 

C 

RLHS*2. •(H1-2.»H2»H3) »(H1-H2-HDELHAX) »«2 
RHS=HDELMAX» (H1-H3) »*2 
IF  (H3.GE.H1.0R.RLHS.GE.RHS)  GO  TO  550 
C 

C  STEP  5 (DEFINE  ZI=A(NMAX) -A<0)  CALCULATE  RLAMBA  SO  THAT 
C  H( A(NMAX) *RLAMBA*ZI )  IS  A  MINIMUM  AND  REPLACE  THE 

C  ZI(IRSAVE)  VALUE  BY  THE  21  VALUE  ABOVE  AND  USE 

C  A(NMAX)*RLAMBA*2I  AS  THE  STARTING  POINT  FOR  THE  NEXT 

C  ITERATION 

C 

IR1  =  1 

391  IF< IRST • EQ • 1 )  THEN 

READ ( 8 , • ) IR1 , HDELM AX , IRSAVE ,<(2I(JfI),J=l. NMAX ) ,1=1, NMAX ) 

IRST=0 

ENDIF 

IF < IR1 .GT.NMAX)  GO  TO  392 
DO  390  IR=IR1,NMAX 

ZI ( IR . IRSA VE ) = A < I R ) - ASAVE ( IR ) 

IF  (ZKIR,  IRSAVE)  .EQ. 0.0)  ZI  (IR.  IRSAVE)  =  .01 
390  CONTINUE 

392  CONTINUE 

RLAMBA=0.0 
DLTA=PERCENT*RLAMBA 
IF  (ABS(DLTA) .LE..05)  DLTA=.05 
RLAMPL=RLAMBA»DLTA 
RLAMMN=RLAMBA - DLTA 
C 

C  ONE  PARAMETER  OPTIMIZATION  FOR  STEP  5  OF  ALGORITHM 
C 

DO  480  K=1 ,KMAX 
RLAMOLD=RL,AMBA 
DO  440  J'l.NMAX 

APLUSDL( J) =A ( J) f RLAMPL»2I (J, IRSAVE) 

AMINDL( J) =A( J) »RLAMMN»2I <J, IRSAVE) 

ABAR( J) =A  < J) ♦ RLAMBA »ZI (J, IRSAVE) 

440  CONTINUE 

CALL  PCAL(APLUSDL, NMAX, PPLUSDL, JMAX) 

CALL  PCAL(AMINDL, NMAX, PMINDL, JMAX) 

CALL  PCAL(ABAR, NMAX, P, JMAX) 

CALL  EXPGRD <  X . Y , P , AA A , CCC , DDD , GGG , WWW , YMIN , 

1  YMAX.IMAX, JMAX, ICOL) 

CALL  EXPGRD  (X, YPLUSDL, PPLUSDL, AAA, CCC, DDD, GGG, WWW, YMIN, YMAX, 

1  IMAX, JMAX, ICOL) 

CALL  EXPGRD  (X. YMINDL, PMINDL. AAA, CCC, DDD, GGG, WWW, YMIN, YMAX, 

1  IMAX, JMAX, ICOL) 

CALL  BLIMP(X,Y,IMAX) JMAX, U, ITERATE, UE.DELTA, RE) 

CALL  BLIMP ( X . YPLUSDL , IMAX , JMAX , UPLUSDL , ITERATE , UE , DELTA . RE) 

CALL  BLIMP(X, YMINDL, IMAX, JMAX, UMINDL, ITERATE, UE, DELTA, RE) 


& 
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CALL  NEWP1 (U , P, UPLUSDL, UMINDL, DUUU , RLAMPL, RLAMHN , 

1  RLAHBA , RMSDUUU ,H,K, PERCENT , IMAX , JMAX , 

2  ABAR2 , NMAX , X , Y , AAA , CCC , DDD , GGG , WWW , YMIN , YMAX , ICOL , 

3  ITERATE , UE , DELTA , P3 , Y3 , D3U , U3 , A , 21 , IRSAVE , WHO , WH2 , WH3) 
WRITE(6,900)  IRSAVE,(A(J> ,J=1,NNAX) 

900  FORMATdX,  'ABOVE  RESULTS  ARE  FOR  OUTER  OPT  AND', 

1  '  IRSAVE*' , 15, / , '  THE  OLD  AU>  VALUES  ARE*', 

2  5E13  5) 

WRITE  (6,901)  *  (ZKJ, IRSAVE)  ,J*1,NMAX) 

901  FORMAT  (IX,'  THE  OLD  ZKJ, IRSAVE)  VALUES*'  ,5E13. 5,/, 

1  IX,'  POLD  IS  CALCULATED  BY  USING', 

2  '  A* A (OLD) ♦ (RLAHBA  OLD)*(ZI  OLD)',/,/) 
ERRORL*ABS(RLAHOLD-RLAMBA) *100.0 

IF  (RLAMOLD.NE.O.O)  ERRORL=ERRCRL/ABS(RLAMOLD>  GO  TO  540 
IF  ( ERRORL . LT . ERMXL )  GO  TO  540 
480  CONTINUE 

540  CONTINUE 

DO  490  J=1 ,NHAX 
A(J)=A(J)+RLAMBA*ZI( J,IRSAVE) 

.  490  CONTINUE 

550  CONTINUE 

370  CONTINUE 

580  CONTINUE 

C 

CC 

CCC 

CALL  POHL ( Y . DELTA , UE , IMAX , JMAX , UEXACT) 

WRITE  (6,400) 

400  FORMAT  (2X, ' I ' ,5X, ' J' , 12X, ' P' , 12X , ' X' , 12X, ' Y' , 12X, 

1  'U',12X, 'UEXACT', 8X,'U-UEXACT') 

DO  610  1=1, IMAX 
IF(I.EG.2.0R.I.EQ.IMAX)  THEN 
DO  613  J=1 , JMAX 

WRITE  (6,500)  I,J,P(J),X(I),Y(I,J),U(I,J), 

1  U£XACT(I,J) ,U(I, J) -UEXACT (I, J)  , 

500  FORMAT  (1X.2I5.6E13.5) 

613  CONTINUE 

ENDIF 

610  CONTINUE 

CALL  PCAL(A,NMAX,P, JMAX) 

CALL  EXPGRD(X,U,P, AAA, CCC, DDD, GGG, WWW, YMIN, 

1  YMAX, IMAX, JMAX, ICOL; 

CALL  BLIMP ( X , Y , IMAX , JMAX , U , ITERATE , UE , DELTA , RE ) 

WRITE  (6,600)  H, RMSDUUU 

600  FORMAT  (IX, 'SUM  OF  THE  SQUARES  OF  DUUU,  H=',E13.5, 

1  /.IX. 'RMS  VALUE  OF  DUUU,  RMSDUUU*' ,E13. 5,/,/) 

CALL  ERROR (U, UEXACT, IMAX, JMAX) 

DO  630  1=2, IMAX 
DO  640  J=1 , JMAX 

ETA=Y( I , J) »SQRT(RE/X ( I ) ) 

WRITE(ll.lO)  ETA,UEXACT(I, J) 

640  CONTINUE 
630  CONTINUE 
10  FORMAT (2E13.5) 

ENDFILE  8 
ENDFILE  9 
ENDFILE  10 
ENDFILE  11 
ENDFILE  12 
ENDFILE  13 
CALL  SECOND  (CPF) 

CPU=CPF-CPI 
WRITE  (6,700)  CPU 
700  FORMAT  (IX, 'CPU*' ,E13.5) 

STOP 


SUBROUTINE  PCAL  (A,NMAX,P,IMAX) 

DIMENSION  A(NMAX> ,P(IMAX) 

DO  10  1=1, IMAX 
C  P(I)=A(1)*I»A(2> 

P( I ) =A(1) »A(2) 

10  CONTINUE 
RETURN 
END 
C 

C  THIS  SUBROUTINE  SOLVES  THE  UNIFIED  DIFF  RELATION  FOR  THE  C  GRID  EQ 
USING  A  TRIDIAGONAL  SOLVER 
C 

SUBROUTINE  EXPGRD(X,Y,P,A,C,D,G,W,YMIN,YMAX,IMAX,JMAX,L) 
DIMENSION  P< JMAX) , A( JMAX) 

DIMENSION  C( JMAX) ,D( JMAX) ,G( JMAX) ,W(JMAX) 

DIMENSION  X(IMAX) , Y(IMAX, JMAX) 

WRITE  (6,100) 

100  FORMAT  (/T15, 'GRID  CAL' ,/,T5, 'USING  A  UNIFIED  DIFFREL') 
Y(L,1)=YHIN 
Y(L, JMAX) =  YMAX 
DO  190  J=2, JMAX-1 
A( J)=EXP( -P( J) ) / <EXP( -P( J) ) +1 .0) 

C(J)=1/(EXP(-P( J) ) ♦ 1 . 0 ) 

D< J) =0.0 
190  CONTINUE 

,G(1)=Y(L,1) 

W(1)=0.0 

DO  290  J=2, JMAX-1 

W(J)=-C(J)/(1.0*A<J)»W<J-1)> 

G(J)=(D(J)*A(J)«G(J-1))/(1.+A(J)*W(J-1)) 

290  CONTINUE 

DC  <90  J=JMAX-1,2,-1 
Y(l,.J)=G(J)-W(J)*Y(L,J+l) 

390  CONTINUE 

DO  490  1=1 , IMAX 
DO  590  J=1 , JMAX 
IF(I.EO.L)  THEN 
Y(I,J)=Y(L„J> 

ELSE 

Y(I,J)=Y(L,J) *SQRT(X(I ) /X(L) ) 

ENDIF 

590  CONTINUE 
490  CONTINUE 

C  WRITE( 14, •)  'ADAPTIVE  GRID  COORDINATES' 

C  WRITE ( 14, 1)((X(I),Y(I,J), 1=1, IMAX), J=l,JNrtX) 

1  F0RMAT(2E13.3) 

WRITE  (6,200) 

200  FORMATdOX, 'TRIDIAGONAL  SUBROUTINE  IS  USED  FOR  GRID  EC') 

RETURN 

END 

C 

C  THIS  SUBROUTINE  CALCULATES  CONVERSION  FOR  PHYSCAL  TO 
C  DORODINITSYN  PLANE.  TRAPEZOIDAL  INTEGRATION. 

C 

SUBROUTINE  TRAPZ ( H , Y , I , J , IMAX , JMAX , SUM , HEND ) 

DIMENSION  H < IMAX.O: JMAX+1) , Y ( IMAX, JMAX) 

SU“-0  0 
DO  10  K-=2, J 
F=Y(I,K)-Y(I,K-1> 

AREA=  F/2» (H(I,K)*H(I,K-1) ) /HEND 
SUM=SUM*AREA 
10  CONTINUE 
RETURN 
END 


no  on 


THIS  SUBROUTINE  CALCULATES  U  VELOCITY  FOR  THE  ITERATION. 


SUBROUTINE  BLIMP<X,Y,IMAX, JMAX, U, IB, UE, DELTA, RE) 

DIMENSION  X ( IMAX) , Y  < IMAX , JMAX ) , U  < IMAX , 0 : JMAX* 1 ) , 

&  V<21 , 19) ,H(21 ,0:60) , W<19) 

DIMENSION  DX(21) ,DY(21,19) ,DYXI(21,19) ,DYB<21,19) 

DIMENSION  ETAX(21 ,19) ,RHO( 19) , RN ( 19) ,T (19) ,DYRE(21 ,19) 

DIMENSION  UNM1(21,19) ,HNM1(21,19) ,VNM1(21,19) ,UIM1(19) 

DIMENSION  HIM1<19) ,UIM2(19) ,HIM2<19) ,Q(21) ,HHREF<21) 

DIMENSION  YP(21,19) ,R(21 , 19) ,DELTA(IMAX> ,DYF(21,19) 

DIMENSION  TAU(21) ,CF<21) ,UD(21 ,19) ,TW<21) 

REAL  MU , MUINF, L,M , MINF ,  MUW . MUEDG , ME, MUSTAR 
DATA  NT , ICNT , KT.NTW , IRST/020 , 01 0 . 10 , 005 , 0/ 

DATA  MINF , PINF , TINF/ .01 , 21 16 . 2 , 530 . / 

DATA  MINF, PINF, TINF/14. 24, .50343,48.6/ 

DATA  DT,L, PR, EPS/200. ,28.95,1. , .00001/ 

DATA  CP, RG, GAM/6006., 1715., 1.4/ 

DATA  FSOR/1 .0/ 

REWIND  13 
REWIND  12 
REWIND  11 
REWIND  10 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  THIS  PROGRAM  SOLVES  THE  UNSTEADY  BOUNDARY  LAYER  EQUATIONS 
C  ALSO  IT  SOLVES  FOR  THE  HEAT  RATE  AT  THE  WALL  AND 
C  DETERMINES  A  HEAT  TRANSFER  COEFFICIENT 
C 


C  SYMBOLS: 

C  AE 

C  AINF 

C  CF 

C  COSJ 

C  CP 

C  DEL 

C  DELM 

C  DELTA 

C  DX 

C  DY 

C  DYXI 

C  DYB 

C  DYC 

C  DYF 

C  DT 

C  EPS 

C  EMAX 

C  GAM 

C  GGM1 

C  H 

C  HC 

C  HE 

C  H.NM1 

C  HIM1 

C  HIM2 

C  HINF 

C  HHREF 

C  HREF 

C  KT 

C  I 

C  ICNT 

C  ICQNP 

C  IMAX 

C  IMET 

C  IPRT 

C  IR 

C  IRST 

C  J 

C  JMAX 

C  JMM1 

C  K 

C  KT 


EDGE  SPEED  OF  SOUND 
INFINITY  SPEED  OF  SOUND 
SKIN  FRICTION  COEFFICIENT 
SOR  CONSTANT 

COEFFICIENT  OF  HEAT(FT-LB/SLUG/R) 

CHANGE  IN  VALUE  AFTER  ITERATION 
MAX  DEL 

BOUNDARY  LAYER  THICKNESS 
DERIVATIVE  OF  X  WRT  XI 
DERIVATIVE  OF  Y  WRT  ETA 
DERIVATIVE  OF  Y  WRT  XI 

Y  DIFFERENCE  BETWEEN  J  AND  J-l 
CENTRAL  Y  DIFFERENCE 

Y  DIFFERENCE  BETWEEN  J*1  AND  J 
DELTA  TIME 

CONVERGENCE  EPSILON 
MAX  ERROR 

RATIO  OF  SPECIFIC  HEATS  (GAMMA) 

GAMMA/GAMMA-1 

TOTAL  ENTHALPY 

HEAT  COEFFICIENT (BTU/FT2/S/R) 

EDGE  ENTHALPY 

TOTAL  ENTHALPY  AT  OLD  TIME  LEVEL 
ENTHALPY  AT  1-1  LOCATION 
ENTHALPY  AT  1-2  LOCATION 
INFINITE  ENTHALPY 
H/HREF  CALCULATED 

REFERENCE  HEAT  C0EFFICTENT(BTU/FT2/S/R) 
THEORETICAL  NONISOTHERMAL  HEAT  COEFF 
COUNTER  IN  XI  DIRECTION 

THE  NUMBER  OF  TIMESTEPS  BETWEEN  PRINTOUTS 

INCOMPRESSIBLE  CASE  FLAG 

MAX  NUMBER  OF  STEPS  IN  XI  DIRECTION 

NUMERICAL  YXI  CALCULATION  FLAG 

FIELD  LISTING  FLAG 

OUTSIDE  DELTA  CALCULATION  FLAG 

RESTART  INDICATOR 

COUNTER  IN  ETA  DIRECTION 

MAX  NUMBER  OF  STEPS  IN  ETA  DIRECTION 

JMAX-1 

COUNTER  IN  ITERATIONS 
MAX  NUMBER  OF  ITERATIONS 


C  L  LENGTH  OF  PLATE  (FEET) 

C  M  AXISYMHETRIC  COEFF. 

C  MINF  MACH  AT  INFINITY 

C  MU  VISCOSITY,  NONDIMENSIQNAL 

C  MUEDG  EDGE  VISCOSITY 

C  MUINF  VISCOSITY  AT  INFINITY(SLUGS/FT/S) 


C  N 

C  NNTS 

C  NT 

C  NTS 

C  NTW 

C  PHI 

C  PINF 

C  PRESS 

C  PRESSO 

C  PR 

C  PE 

C  0 

C  RE 

C  REE 

C  RG 

C  RHO 

C  RHCE 

C  RM 

C  RE 

C  RINF 

C  T 

C  TAU 

C  TEO 

C  THETA  ■ 

C  TIME 

C  TE 

C  TINF 

C  TOINF 

C  TS 

C  TW 

C  TW1 

C  TW2 

C  U 

C  UNM1 

C  UE 

C  UIM1 

C  UIM2 

C  UINF 

C  V 

C  VNM1 

C  W 

C  XTW2 


COUNTER  IN  TIME 

MAX  NUMBER  OF  TIME  STEPS 

NUMBER  OF  TIME  STEPS  TO  CALCULATE  OMEGA 
NONISOTHERMAL  CORRECTION  FACTOR 
PRESSURE  INFINITY 
NONDIMENSION ALI2ED  EDGE  PRESSURE 
NONDIMENSIONALIZED  INITIAL  PRESS 
PRANDTL  NUMBER 

PRESSURE  AT  THE  EDGE  OF  B.L.U.B/FT2) 

HEAT  RATE  AT  THE  WALL(BTU/FT2/S) 
REYNOLD'S  NUMBER  ON  L 
EDGE  REYNOLDSS  NUMBER  ON  L 
GAS  CONSTANT  (FT-LB/SLUG-R) 

DENSITY 

EDGE  DENSITY 

DENSITY  TIMES  VISCOSITY 

REYNOLD'S  NUMBER 

DENSITY  AT  INFINITY  (SLUG/FT»»3) 

TEMPERATURE 

SKIN  FRICTI0N(LB/FT2) 

EQUIVALENT  TEMPERATURE(R) 

WEDGE/CONE  HALF  ANGLE 
TOTAL  TIME  ACCUMULATED 
TEMPERATURE  AT  THE  EDGE  OF  B.L.(R) 
TEMPERATURE  AT  INFINITY  (DEG  RANKINE) 
TOTAL  TEMPERATURE  AT  INFINITY  (RANKINE) 

TEMPERATURE  AT  THE  WALL(R) 

TEMPERATURE  AT  WALL  BEFORE  PANEL 

TEMPERATURE  AT  WALL  AT  PANEL 

VELOCITY  IN  X  DIRECTION 

VELOCITY,  X  DIRECTION  AT  OLD  TIME  LEVEL 

VELOCITY  AT  THE  EDGE  OF  B.L.(FT/S> 

VELOCITY  VECTOR  AT  1-1  LOCATION 

VELOCITY  VECTOR  AT  1-2  LOCATION 

VELOCITY  AT  INFINITY  CONDITION  (FT/SEC) 

VELOCITY  IN  Y  DIRECTION 

VELOCITY, Y  DIRECTION  AT  OLD  TIME  LEVEL 

OPTIMIZATION  OMEGA 

POSITION  OF  PANEL  ON  WEDGE 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 


C  WRITE(6, * ) ' NUMERICAL  YXI  METRIC  AFTER  IMET= ' 
READ(5, » ) IMET 

C  WRITE (6,*) ' L (FT)  =  ,  PR=' 

READ ( 5,  » )  L ,  PR 

C  WRITE (6, * )  'INPUT:  CONE=l . ,WEDGE=0. ' 

READ(5, • )  M 

C  WRITE(6,»>  ' IR=0,R=F (X)  ONLY' 

C  WRITE(6, • )  ' IR= ' 

READ (5, *)  IR 

C  WRITE ( 6, *> 'MINF.PINF, TINF*' 

READ(5,»> MINF, PINF, TINF 

C  WRITE (6, * ) ' TWI , TW2, XTW2* ' 

READ(5.*)TW1,TW2,XTW2 

C  WRITE(6, * ) ' IPRINT=0  NO  LISTING  ON  FILE  FIELD' 

C  WRITE(S, * ) ' IPRT* ' 

'  READ(5, *) IPRT 

C  WRITE (6 , * ) 'NT, ICNT.KT, IRST,NTW= ' 

READ(5, »)NT, ICNT,KT,IRST,NTW 

C  WRITE (6, • ) 'DT,EPS=' 

READ(5, « )DT,EPS 


C  WRITE<6.*)  'THETA, BETA  GUESS?' 

READ (5,*)  THETA, PO 
PI=AC0S<-1. ) 

JMMWMAX-1. 

GGM1=GAM/ (GAM-1.) 

HUINF*(TINF««1 . 5»2 .2685E-8) / <T1NF*198 .6) 

AINF  =SQRT<GAM»RG»TINF) 

UINF=MINF«AINF 
RINF=PINF / (RG»TINF) 

TOINF=TINF»  < 1* (GAM-1 . )«.5»MINF»«2> 

HINF  =CP»TINF*UINF«»2/2 
COSJ=COS ( PI/FLOAT ( JMAX > ) »FSOR 
RE=RINF*UINF*L/MUINF 
THETA^THETA'PI/iaO. 

IF (THETA.NE.O.O.AND.MINF.GE.l .0)THEN 
CALL  BETA ( THETA , MINF , P I NF , TI NF , PO , CP , RG , ME , PE , TE , 

&  HE, UE) 

A£=SQRT(GAM»RG«TE) 

HUEDG  =  (TE»»1.5«2. 2S85E - 8 ) / ( TE* 198 . 6 ) 

RHOE=PE/(RG«TE) 

REE=RHOE*UE*L/MUEDG 

ELSE 

A£=AINF 

HE-MINF 

PE=PINF 

TE=TINF 

HE=HINF 

UE=UINF 

RHOE=RINF 

REE=RE 

MUEDG=MUINF 

ENDIF 

HEND=HE/UINF**2 
IF( IB.EQ.O)  THEM 
WRITE<9, • ) 'RE= ,PR=  ',RE,PR 

WRITE (9,  *)  'RINF,UINF,MUINF='  ,RINF,UINF,MUINF 
WRITE (9, • ) 'MINF , AINF ,TINF, PI NF  »' .MINF , AINF ,TINF,PINF 
WRITE (9,*) 'ME,PE,TE, AE,UE=' ,ME,PE,TE, AE.UE 
WRITE(9, * )  ' NT, ICNT,KT,IRST,DT,EPS,NTW , IMET= , ' , 

1  NT , ICNT , KT , IRST , DT , EPS , NTW , IMET 

ENDIF 

PRESS=PE/ (RINF»UINF**2) 

PRESSQ=PRESS 

HREF=-332»CP*SQRT(NUINF«RINF«UINF)/PR»*(2./3.) 

C 

C  GENERATE  A  GRID  AND  INITIAL  CONDITIONS 
C 

CALL  FIRST<IMAX, JMAX,REE,CP,TW,TE,UE,UINF,X, Y,U,V,H,PR, 
8.  XTW2.TW1.TW2, IMET. PE, RG.M.RHOE, MINF, DELTA, L) 

IF  (Y<IHAX, JMAX) . LT. DELTA (1MAX) , AND.MINF.GT. .05) 

1  GO  TO  920 
DO  15  1  =  1  ,IMAX 
Ud,0)=0. 

H< I ,0) =0.0 
Ud,JMAX*l)=UE/UINF 
15  • Hd ,  JMAX«1 ) =HEND 

C 

DO  20  1=1 , IMAX 
U< I. JMAX) =UE/UINF 
20  H< I , JMAX) =HEND 

C 

C  CALCULATE  METRICS 

CALL  METRIC  < X , Y , RE , I MET , DX , D Y , DYXI , DYB , DYC , DYF , ETAX . 

&  DYRE, IMAX, JMAX, JMM1) 

C 

C  PUT  INFORMATION  INTO  THE  OLD  TIME  LEVEL  INITIALLY 
C 

DO  30  1=1, IMAX 
DO  30  J=1 , JMAX 

UNM1 ( I , J)=U< I , J) 
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VNM1(I,J)=V(I, J) 

30  HNM1<I,J)=H<I,J> 

DO  40  J=1 , JMAX 
IF  (J.LE.ll)  THEN 
W<J)=1.6-<J-1)*.06 
ELSE 
W<J>=1.0 
ENDIF 

40  CONTINUE 

C 

C  START  THE  TIME  LOOP 

TIME=0.0 
DO  100  N=1 ,NT 
■  TIME*TIME*DT 

DPDT  =  <  PRESS-PRESSO ) 

PRESSO=PRESS 

C 

C  START  BACK  AT  BEGINNING  OF  GRID- -REINITIALIZE  VECTORS 

C 

120  DO  130  J*1 , JMAX 

UIM1 < J) =0.0 
130  HIM1(J>=0.0 

C 

CALCULATE  THE  OPTIMIZATION  FACTOR  OMEGA  EVERY  NTW  TIME  STEPS 
C 

IF(N.EO. (N/NTW>  »NTW>THEN 

CALL  MURHO  <  H , U , GGM1 , PRESS , UINF , CP , MUINF , IMAX , JMAX ,RM . 

1  W.TE) 

CALL  OMEGA  <  U , V , DT , DX , DY ,  DYB ,  DYF , DYRE , ETAX , RM , W , WN , 

&  IMAX, JMM1, JMAX, COS J) 

ENDIF 

IF(IR.EO.l)  CALL  RADIUS<X, Y, DELTA, IMAX, JMAX, THETA, R,IR> 
C 

C  START  I  LOOP- -MARCH  IN  XI  DIRECTION 
C 

DO  200  1=2. IMAX 
IM1=I-1 
IM2=I-2 
DELM=0. 0 

C  PUT  NEW  INFORMATION  INTO  THE  VELOCITY  AND  ENTHALPY 
C  VECTORS  AT  THE  1-1  AND  1-2  LOCATIONS 
C 

DO  210  J=1 , JMAX 
UIM2(J)=UIM1(J) 

UIM1 ( J) =U<I-1 , J) 

HIM2( J) =HIM1 ( J) 

210  HIM1(J)=H(I-1,J) 

C 

C  DEFINE  CONSTANTS  TO  DETERMINE  THE  COEFFICIENTS  FOR  THE 
C  U  AND  V  VELOCITY  EQUATIONS  AND  THE  ENTHALPY  EQUATION 
C 

IF(I.EQ.2)THEN 

CIJ1=1. 

CIM1=1. 

CIM2=0. 

CVI J=- .5 
CVIM1=.5 
. CVIM2=0. 

ELSE 

CIJ1=1 .5 
CIM1=2.0 
CIM2=-.5 
CVIJ=-.75 
CVIM1=1 , 

CVIM2=-.25 

ENDIF 

C 

C  START  THE  ITERATION  LOOP 
C 

DO  300  K=1 ,KT 


•v  .-'.--j. 
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EMAX=0.0 

EUMAX=0. 

EVMAX=0. 

EHMAX=0. 

CALCULATE  NEW  VALUES  FOR  TEMPERATURE, DENSITY,  AND  VISCOSITY  USING 
THE  NEW  VALUES  OF  U,  V,  AND  H  (CALCATED  AT  THE  OLD  ITERATION  LEVEL) 

DO  310  J*1,JNAX 

T(J)=H(I,J)-U(I,J)»*2/2. 

TEMP=T< J) »UINF»«2 . /CP 
RH0<J)=GGM1*PRESS/T<J) 

IF (TE.LT.200. )  THEN 

MU=TEMP»» .5*2.32E-8/ ( 1 . *220./ (TEMP»10. »• <9.^TEMP> ) ) /MUINF 
ELSE 

HU=  <TENP*«1 .5»2.2685E-8> / (TEMP»198 .6) /MUINF 
ENDIF 

RM< J)*RHO< J) »MU 

START  J  LOOP--MARCHING  IN  ETA  DIRECTION 

DO  400  J=JMM1 ,2,-1 
JM1= J-l 
JM2= J-2 
JP=J*1 
UOLD=U<I,J> 

HOLD-H(I, J) 

IF(K.EQ.1)THEN 
W0PT=1.0  . 

ELSE 

WOPT*W(J> 

ENDIF 

SETUP  CONSTANTS  TO  DETERMINE  COEFFICIENTS  IN  U,  V,  AND  H  EQUATIONS 

SOLVE  FOR  U.V.AND  H 

DRMF= (RM< JP) »RM ( J) ) /DYF ( I , J) » .5 
DRMB=  (RM( J) ♦RM( JM1 ) > /DYB(I , J) • . 5 
CON=ETAX<I,J)*U(I,J)  .♦  V(I,J)/DY(I,J) 

IF ( CON. LT. 0.0) THEN 
IF(J.EQ.JMNl)  THEN 
CIJ2=-1.0 
CJM1  =  -1 .0 
CJM2=0.0 
ELSE 

CIJ2=-1 .5 
CJMl=-2.0 
CJM2=  .5 
ENDIF 
DRH1=DRHF 
DRM2=DRMB 
UWP 
J2-J*2 
J3= JN1 
ELSE 

IF(J.EQ.2)THEN 
Cl J2=l . 

CJM1=1. 

CJM2=0. 

ELSE 

CIJ2=1 .5 
CJM1=2.0 
CJM2=- .5 
ENDIF 
DRM1=DRMB 
DRM2=DRMF 
J1=JM1 
32-JH2 
J3=JP 
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ENDIF 

C  SOLVE  FOR  U  VELOCITY 

CIJ=l.  *  CIJ1»DT/DX(I)»U<I,J>  ♦  Cl J2»DT»C0N 
&  ♦  DT»DYRE(I , J) • <DRMF*DRHB) 

C0NJ1=CJM1»DT»C0N  *  DT»DYRE< I , J> »DRM1 
CON J2=CJM2»DT  »CON 
CONIN=DT/DX<I)*U(I,J> 

C0NVS=DT»DYRE<I,J)«DRN2 

US= (UNMl ( I , J)  ♦  CIM1»C0NIM«UIH1 < J)  ♦  CIM2«C0NIM*UIM2< J>  ♦ 

&  C0NJ1»U(I ,vil)  *  C0NJ2»U<I,J2>  ♦  CONVS»U<I . J3> > /CIJ 
U(I,J)=WOPT»US  *  <1. -WQPT) •U(I , J) 

C  SOLVE  FOR  ENTHALPY 

CON=ETAX(I , J) «U ( I , J)  ♦  V(I,J)/DY<I,J) 

IF ( CON. LT. 0.0) THEN 
IF(J.EQ.JNNl)  THEN 
CIJ2=-1.0 
CJN1*-1 .0 
CJM2=0.0 
ELSE 

CIJ2*-1.5 
CJNl=-2.0 
CJN2* .5 
ENDIF 
DRH1=DRNF 
DRN23DRMB 
J1=JP 
J2=J*2 
J3=JH1 
ELSE 

IF<J.EQ.2)THEN 
CIJ2=1 i 
CJN1=1 . 

CJH2=0. 

ELSE 

0112*1.5 
CJN1=2.0 
CJM2=-  .5 
ENDIF 
DRM1=DRHB 
DRM2=DRNF 
J1=JX1 
J2= JN2 
J3=JP 
ENDIF 

CIJ=1.  *  CIJ1*DT/DX(I)»U(I,J)  ♦  Cl J2»DT»C0N 
S.  *  DT»DYRE< I , J) » <DRMF*DRKB) /PR 

C0NJ1=CJM1»DT*C0N  »  DT*DYRE( I , J) •DRN1/PR 

C0NJ2=CJM2*DT»CQN 

CONIM=DT/DX  < I ) «U ( I , J) 

CONVS=DT»DYRE(I , J) «DRN2/PR 
TCON=DYRE( I , J) *DT  »<(PR-1.)/PR>*,5 

HS=<HNM1(I,J)  ♦  Cl'Ml*CONIM*HIM.l  ( J>  ♦  CIM2»C0NIM»HIM2( J)  ♦ 

&  C0NJ1»H<  I ,  J1 )  ♦  C0Mo'2»H ( I ,  J2)  ♦  C0NVS»H(I,J3)  * 

&  DPDT/RHO(  J)  *TC0N*  (DRMF  «U< I , JP)  »*2- (DRMF  *DRMB)  *U(I,J)*»2* 

S.  DRMB»U<I,JM1)»«2)>/CIJ 

H(I,J)=«OPT*HS  ♦  < 1 . -WOPT) *H( I , J  > 

COMPUTE  THE  ERROR 

EU=ABS(U< I , J) -UOLD) 

£H=ABS(H(I, J)-HQLD) 

IF(EU.GE.EUMAX)EUMAX=EU 
IF  (EH.GE.EHMAX)EHMAX=EH 
IF(£UMAX.GE.£MAX)EMAX=EUHAX 
IF(EHMAX. GE.EKAX) EMAX=EHNAX 

END  J  LOOP 
400  CONTINUE 
C  CALCULATE  V  VELOCITY 
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DO  320  J*2, JMM1 
JMW-1 
VOLD=V(I, J) 

C  IF  (IR.EQ.O)  THEN 

V<I,J)*V(I,JM1)*<.5» (DYXI (I,J)+DYXI(I,JM1)>* 

1  <U(I,J>-U<I,JM1))+<Y<I,J>-Y(I,JM1)>* 

2  <CVIJ*(U(I,JKU<I,JM1>)* 

3  CVIM1»<UIM1<J)*UIM1<JM1)K 

4  CVIM2«  <UIM2( J) +UIN2  < JN1 ) ) ) ) /DX ( I ) 

C  ELSE 

C  V(I,J)=<1/R<I,J)»«N)»(R(I,JM1)»«M»V<I,JM1)  ♦  R<I,J)»*M 

C  &»VCNI J*U(I , J)*R(IM1, J>  *«M*CVIM1»D1*UIM1 < J)  ♦  R(IM2,J)**M 

C  &»CVIM2»D1»UIM2(  JKRd ,  JM1 ) •*N«VCNJ1»U ( I , JM1) *R( IM1 , JM1) «*M 
C  S.»CVIM1*D1»UIM1 ( JM1) *R< IM2, JM1) »*M»CVIM2»D1*UIM2( JM1 ) ) 

C  ENDIF 

EV= ABS ( V ( I , J ) -VOLD  > 

IF<EV.GE.EVMAX>EVMAX=EV 
IF(EVMAX.G£.EMAX>£MAX=EVMAX 
320  CONTINUE 

C 

C  CHECK  CONVERGENCE 

C 

IF (EMAX.LE.EPS)  GO  TO  420 
C 

C  END  ITERATION  LOOP 

300  CONTINUE 

420  DO  430  J=2, JMM1 

DELrABS(U(I , J)  -  UNM1(I, J)  ) 

IF(DEL.GT.DELM)DELM=DEL 
DEL=ABS(H<I, J)  -  HNMKI,  J) ) 

I F < DEL . GT . DELM ) DELH > DEL 
430  CONTINUE 

V  ( I',  JMAX )  = V  ( I ,  JMN1 ) 

DO  330  J=1,JMAX 
U(1,J)=U(2,J) 

H(1.J)=H<2,J) 

330  V(1.J)=V(2.J) 

C  MOVE  THE  RESULTS  INTO  THE  OLD  TIME  LEVEL 
C 

DO  440  J=1.JMAX 
UNMKI,  J)=U(I,  J) 

HNMKI. J)=H<I,J) 

440  VNMKI,  J)=V(I,  J) 

C 

C  CALCULATE  THE  HEAT  RATE  AT  THE  WALL 
C 

FDEL Y  *(-3.»Y(I,l)  ♦  4.*Y(I,2)  -  Y(I,3)>/2. 

IF  ^TE.LT.200. )  THEN 

MUW=2.32E-8»TW< I) ** .5/(1. *220,/ <TW(I) »10. •• <9. /TW< I ) ) ) ) 
ELSE 

MUW=<TW<I)#«1.5*2.2685E-8>/CTW<I)*198.6) 

ENDIF 

C  CONDUCTIVITY  IS  EQUAL  TO  MUW»CP/PR=K 

Q( I) =MUW»RHO( 1) »(-3.»T(l)  ♦  4.«T(2)  -  T(3))«  RINF 
S.  »UINF»*2/(2. »L»PR*FDELY»RHOE> 

C 

C  TAW  SHOULD  USUALLY  BE  BASED  ON  EDGE  CONDITIONS  FORMALLY 
C  BUT  CAN  BE  BASED  ON  FREESTREAH  IF  YOU  ARE  CONSISTENT!!!! 

C 

TAW=TE» ( 1 . ♦SORT (PR) » (GAM-1 . ) • .5«ME»*2) 

HC=Q( I ) / (TAW  -  TW < I } ) 

HHREF ( I >  =HC/ HREF  »SQRT ( X ( I >  «L  > 

TAU(I) =MUW*UINF« ( -3»U( I , 1) *4»U(I,2) -U(I ,3) ) • .5 
1  /(L*FDELY»RH0(1)) 

CF(I ) *2. »TAU (I) / (RINF«UINF»»2) 

C 

C  END  I  LOOP 
200  CONTINUE 
C 

C  WRITE  OUT  SOLUTION  AT  A  TIME  STEP 
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c 

IF(N.EQ.l .OR.N.EQ. <N/ICNT)»ICNT)THEN 
WRITE(10,1)N,TINE 

WRITE ( 6 , • ) ' DELT  = ' , DELM ,'K=',K,'U,V,H  ERR®' , EUMAX . EVMAX , EHMAX 
WRITE < 10, •) 'DELT®' , DELM, 'K*',K,'U,V,H  ERR®' .EUMAX, EVMAX, EHMAX 
WRITEO,  12) 

DO  150  I®1 , IMAX 

IF  < I . £Q . 2 . OR . I . EQ . IMAX )  THEN 

WRITE (10, 2) I ,X(I) 

IF(IPRT.E0.0)G0  TO  181 
WRITE(10,3> 

DO  160  J*1 , JMAX 
T(J)*H(I,J)-U<I,J)»®2/2. 

TEMP*T< J)»UINF»»2/CP 

WRITE<10,4) J, Y(I, J) ,U(I,J),H(I,J),T(J) ,TEMP,V(I , J) ,RHO<J) 

160  CONTINUE 
ENDIF 

161  IF(I.EQ.l)  GO  TO  150 
WRITE<10,5)  Qd)  ,HHREF(I> 

TIHEP=TIME*L/UIMF 

C  CF  RATIO  TO  BLASIUS  CAN  BE  BASED  ON  EDGE  OR  FREESTREAM  RE 
CFCF=CF(I)/.664»SQRT<REE«X(I) ) 

WRIT£(9,7>TIMEP,Xd)  ,CFCF,HHREF(I),TW<I) 

IF(MINF.LE. .05. AND.N.EQ.NT)  THEN 
WRITE<13,10)  X(I).CFCF 
WRITE(12,10)  X(I),HHREF<I) 

ENDIF 

150  CONTINUE 
ENDIF 

C  END  TIME  LOOP 

100  CONTINUE 
C 

C  COMPUTE  PHYSICAL  Y  VALUES 

IF  (MINF.LE- .05)  THEN 
DO  700  1=2. IMAX 
DO  710  J=1 , JMAX 
ETA= Y (I , J) -SORT ( RE/X (I ) ) 

WRITE* 11. 10)  ETA.U(I,J) 

710  CONTINUE 
700  CONTINUE 
GO  TO  920 
ENDIF 

DO  BOO  J=1 , JMAX 
YP(1, J)=0.0 

C  WRITE <11, 10)X(1),YP<1,J) 

800  CONTINUE 
C  WRITli<13.11) 

DO  810  1=2, IMAX 
DO  820  J=1 , JMAX 
T( j ) =H(I, J)  -  U(I,J>**2/2 
820  RHO( J) =GGM1 *PRESS/T ( J) 

C 

YP(I,1)=0.0 

C  WRI.TE <  11, 10)Xd)  ,YP(I,1) 

DO  b30  J=2, JMAX 

C  YP<I,J)*YP(I,J-1)  ♦  (1/RH0(J)  *  1 /RhO< J-l ) )/2»<Y(I,J)-Yd,J-l)) 

CALL  TRAPZ(H, Y, I, J, IMAX, JMAX, SUM, HEND) 

YPCI , J) =SUM 

ETAC® YP (I , J ) .SORT  <  UE«RHOE*L/ ( MUEDG»X ( I) ) ) 

C  WRITE  (11 ,10)Xd),YP(I,J) 

T(J)=T(J)»UINF**2./(CP*TE) 

UD(I,J)=U(I,J)®UINF/UE 
WRITE(11,10)  ETAC,T( J) 

WRITE(13,10)  ETAC.UDd, J) 

C  WRITE (13,8)  T(J) ,U(I, J) ,ETAC 

830  CONTINUE 
810  CONTINUE 

DO  825  1=1, IMAX 
DO  826  J=1 , JMAX 
WRITE(11,10)  X(I) , Y(I, J) 
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826  CONTINUE 
825  CONTINUE 
C 

C  THIS  SECTION  FIGURES  OUT  THE  THEORETICAL  NONISOTHERMAL  VALUES 
C  FOR  H/HREF 
C 

C  CAUTION!  TSTAR  BASED  ON  TW1!  DO  YOU  WANT  CONSTANT  REFERENCE  ■ 

C  OR  ACTUAL  TW? 

840  TAW*TE*(1. ♦SORT (PR) • (GAN-l. ) * .5*ME**2) 

TSTAR=TE* .5* (TW1-TE) ♦ .22* (TAW-TE) 

RSTAR*PE/ (RG*TSTAR) 

IF (TE.LT.200. >  THEN 

MUSTAR* ( 2 . 32E-8»TSTAR»* .5/(1. *220 . / <  TSTAR* 10 . • • (9 . /TSTAR ) ) ) ) 

ELSE 

MUSTAR*  (TSTAR*  *1. 5*2. 2685E-8)/ (TSTAR*  198. 6) 

END  IF 

DO  900  I*1,IMAX 

HIS0*.332»CP*SGRT(MUSTAR«RSTAR*UE)*SQRT(X(I)*L)/(PR**(2./3.)*HREF) 
WRITE(9, •)  X(I) , 'HISO*' , HISO 
900  CONTINUE 

920  CONTINUE 

1  FORMAT (///,2X,'N*' ,I6,5X,'TIME*',F15.6,1X, 'SEC' ) 

2  F0RNAT(//,2X, 'I*' ,13, SX, 'X-LOCATION*' .F13.5) 

3  F0RMAT(/,3X,'J',7X,'Y-L0C' ,10X,'U'  ,11X, 'H'  ,11X, 'T' ,11X, 
fc'TEMP(R) '  ,  11X,  'V',HX,  'RHO' ) 

4  FORMATdX,  I3.5F13.5.2F15.8) 

5  F0RMAT(2X, 'HEAT  RATE  (Q>  =  ' , F13. 6, 5X, 'H/HREF*' ,F13.6) 

6  FORMAT ( 'THIS  SOLUTION  CONVERGED  WITH  KT  ITERATIONS') 

7  F0RMAT(2X,F15.5,4(2X,F13.5) ) 

8  FORMAT ( 3 ( 2X , £13 . 5 ) ) 

9  F0RMAT(34X,2(2X,F13.5) ) 

10  FORMAT ( 2E13 . S ) 

11  F0RMAT(5X, 'T/TINF' ,9X. 'U/UINF' ,11X, 'ETAC' > 

12  FORMAT ( / . 7X . ' TIME ',15X , ' X ' .  13X , ' CF/CF ' , UX . ' HHREF ' , 

1  10X, 'TWALL' ) 

REWIND  S 
RETURN 
END 
C 

C  THIS  SUBROUTINE  CALCULATES  BETA  AND  EDGE  CONDITIONS 

C 

SUBROUTINE  BETA (THETA , NINF , PINF , TINF , PO , CP , RG , ME , PE , TE , 

&HE.UE) 

C  CALCULATES  EDGE  CONDITIONS  GIVEN  DEFLECTION  ANGLE 
C  SHOCK  ANGLE  CALCULATED  TO  1.0E-08  TOLERANCE 
R£AL*  NINF  NE  HE 

Q(X) 3 ( (GAM*1 ! )/2. *SIN(X) *SIN (THETA) /COS (X-THETA) ♦ 

6(1 . /MINF**2. ) ) **.5-SIN(X) *X 
PI=AC0S(-1.) 

GAM*1.4 
T0L=1 .0E-08 
P0=P0*PI/180. 

Y=Q(P0) -PO 
DO  10  1=1,50 
XNEW=Q(PO) 

Y=Q(XNEW) -XNEW 
ER0R=ABS(XNEW-?0) /ABS(XNEW) 

IF ( EROR . LT . TOL . OR . ABS ( Y ) . LT . TOL)  GO  TO  30 
PO=XNEW 
10  CONTINUE 

30  ME=((2. ♦(GAM-1. >*(HINF*SIN(XNEW))**2.)/((2.*GAH*(HINF* 

&SIN(XNEW) ) **2. - (GAM-1 . ) ) • (SIN(XNEW-THETA) )*»2.))**.S 
TE=TINF* (2. ♦ (GAM-1 . ) • (MINF* SIN (XNEW) ) **2. ) / (2. 

&  *(GAM-1 . ) *(ME*SIN (XNEW -THETA) ) **2. ) 

UE=ME*SQRT(GAM»RG»TE) 

HE=CP»TE*UE**2/2 

PE* (1* (2*GAM/ (GAM*1 ) ) » ( (NINF*SIN (XNEW) )**2-l))*PINF 

RETURN 

END 
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C  THIS  SUBROUTINE  TAKES  THE  BLASIUS  SOLUTION,  MODIFIES  IT  BY 
C  THE  DENSITY,  AND  FORMS  A  GRID.  IT  ALSO  DETERMINES  INITIAL 
C  CONDITIONS 

SUBROUTINE  FIRST(IHAX,JHAX,RE,CP,TW,TE,UE.UINF. 

S,  X,Y,U,V,H,PR, XTW2 , TW1 , TW2 , IMET , PE , RG , M , RHOE , MINF , DELTA , L ) 
DIMENSION  X<IMAX) , Y(IMAX, JMAX) , U(IMAX,0: JHAX+1) ,HdMAX,0:  JMAX+1) 
DIMENSION  VdMAX,  JMAX> , TW(IMAX) , DELTA (IMAX) 

REAL  H,HINF,MUW,L 
C 

C  TO  USE  ANALTICAL  METRIC  FOR  YXI  FOR  1=1  TO  IMET 
C  Y  POINTS  AT  I*IMET  ARE  SCALED  BY  SQRT(X) 

C  THIS  WILL  NOT  GUARANTEE  A  SMOOTH  GRID  PAST  I=IMET 
DO  30  1=1, IMET 
DO  30  J=1 , JMAX 

Yd,J)=Y< IMET,  J)  »SQRT(Xd)  /X ( IMET) > 

30  CONTINUE 
C 

C  DEFINE  THE  WALL  TEMPERATURES--  IF  TW1=TW2  NO  TEMP  STEP 
C  IF  TW1<TW2  TEMP  STEP 

C 

DO  40  1=1, IMAX 
IF(X(I) .LE.XTW2)THEN 
TW(I)=TW1 
ELSE 

TW(I)=TW2 

ENDIF 

40  CONTINUE 

HE=<CP«TE*UE»»2./2. )/UINF»*2. 

DELTA <1) =0.0 
DO  50  I =2, IMAX 
HW=CP=TW  <I)/UINF»»2 
HSUB=HE  -  HW 
IF  (MINF.LE. .05)THEN 
DELTA d ) =5.2/SQRT(RE) »SQRT(X < I ) ) 

DELTAT= DELTA  < I ) / ( 1 . 026*PR • * ( 1 . / 3 . ) ) 

ELSE 

RHOW=PE/(RG*TW(I)) 

IF (TE.LT.200. )  THEN 

MUW=2.32E-fl*TWd)  =  ».5/d .♦220./(TU(I)»10.«#(9./TW(I) > ) ) 

ELSE 

MUW=  <  TW  d  )  •  •  1 .  5= 2 . 2685E -8)/(TW(I)*198.6> 

ENDIF 

DELTA ( I ) * (280»RHOW»MUW«X<I)/ ( (13. •RHOE«RHOE»L*UE» (1»2»M))))»*.5 
DELT  AT = DELT  A ( I) 

ENDIF 

DO  60  J=l, JMAX 
IF(Y<I,J) . LE . DELTA ( I ) ) THEN 

Ud,J)=UE/UINF»d.5*Yd,J>/DELTAd>  -  .5*  ( Yd ,  J) /DELTA ( I) ) »*3. ) 
ELSE 

Ud ,  J)  =UE/UINF 
ENDIF 

V<I,J)=0.0 

IF(YCI.J) .LE.DELTAT)THEN 

H(I,J)=HW  ♦  HSUB* (1.5=Y(I,J) /DELTAT  -  .5= (Y ( I , J) /DELTAT) **3) 

S.  ♦  <U(I,J)/UINF>»»2./2.  *  (UE/UINF) »»2. /2» 
t  ( 1 .5»Y( I , J) /DELTAT- .5* < Y < I ,J) /DELTAT) » *3. ) 

ELSE 

H(I,J)=HE 

ENDIF 

60  CONTINUE 
50  CONTINUE 

DO  70  J=1 , JMAX 
U(l, J)=U<2, J) 

V(1,J)=0.0 
H(l, J)=H<2, J) 

70  CONTINUE 
RETURN 
END 
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COMPUTE  THE  DENSITY  TIMES  VISCOSITY  FOR  EVERY  J  LOCATION 


SUBROUTINE  MURHO (H , U.GGM1 , PRESS , UINF , CP , MUINF , IMAX , JMAX , 

&RN  W  TE) 

DIMENSION  H ( IMAX , 0 : JMAX ♦ 1 ) ,U < IMAX , 0 : JMAX ♦ 1 ) , RM < JM AX ) , W ( JMAX ) 
REAL  MU, MUINF 
DO  10  1*2, IMAX 
DO  20  J*l, JMAX 
TEMP*H(I,J)-U(I,J)**2/2. 

DEN=GGM1»PRESS/TEMP 
T£MP*TEMP»UINF*»2/CP 
IF(TE • LT. 200 • ) THEM 

MU* (2.32E-8»TEMP»* .5/ < 1 .+220. / (TEMP* 10. *• <9. /TEMP) ) ) ) /MUINF 
ELSE 

MU* (TEMP** 1 .5»2.2685E-8) /(TEMP+198. 6) /MUINF 
ENDIF 

RM( J) * DEM* MU 
20  W ( J ) * 1 . 85 

10  CONTINUE 
RETURN 
END 


C  THIS  SUBROUTINE  CALCULATES  THE  SOR  OMEGA 
C 

SUBROUTINE  OMEGA ( U , V , DT , DX , DY , DYB , D YF , DYRE , ETAX , RM , 

&W , WN , IMAX , JMM 1 , JMAX , COS J ) 

DIMENSION  U(IMAX,0:JMAX+1),V(IHAX, JMAX), DX(IMAX) ,DY(IMAX, JMAX) 
DIMENSION  DYB < IMAX , JMAX ) , DYF ( IMAX , JMAX ) , DYRE ( IMAX , JMAX) 
DIMENSION  RM( JMAX) ,W( JMAX) , ETAX (IMAX, JMAX) 

REAL  LAM 
DO  10  1*2. IMAX 
DO  20  J*2,JMM1 
JP=J*1 
JM*J-1 

DIV=1./DT  ♦  1.5«<U(I,J)/DX(I)  ♦  ABS( U (I , J)*ETAX ( I , J)  ♦' 

&  V<I,J)/DY<I,J)>)  ♦  DYRE(I,J)*.5* 

t  (<RM<JP)-RM(J))/DYF(I,J)  *  ( (RM(J)*RM(JM) )/DYB(I,J) ) ) 

LAM  = ( - 2 . »SQRT (ABS(ABS<(ETAX(I,J)»U(I,J)  ♦  V(I , J) /DY(I, J) )  ♦ 
fc  .5*DYRE(I , J) *<  <RM( J) ♦RM ( JM) )/DYB (I ,J) ) )  « 

&  < .5»DYRE(I , J) »<(RM(JP)*RM(J) )/DYF(I, J) ) ) ) ) ) /DIV*COSJ 

IF ( LAM . GT . 0 . 97 ) LAM=0 . 97 
WN=2./(1.*SQRT(1.-LAM»*2) ) 

IF(WN.LT.tf(J) )W( J) =WN 
20  CONTINUE 

10  CONTINUE 

RETURN 
END 
C 

C  COMPUTES  H3  FOR  A  CONE- ASSUME  CONSTANT  TRANSVERSE  RADIUS 
C  ASSUME  YC05< ALPHA)  SMALL.  IF  NOT  DELETE  NEXT  LINE. 

C 

SUBROUTINE  RADIUS(X,Y, DELTA, IMAX, JMAX, THETA, R,IR) 

DIMENSION  X(IMAX),Y(IMAX, JMAX), R(IMAX, JMAX) .DELTA (IMAX) 

DO  10  1=1, IMAX 
DO  20  J=1 , JMAX 
IF( Y(I , J) . LE. DELTA ( I ) )  THEN 
R(I,J)=  X(I) »SIN(THETA) 

ELSE 

IF(IR.EO.O)  THEN 
R(I, J)=X(I) »SIN(THETA) 

ELSE 

R(I,J)=X(I) »SIN (THETA ) *Y<I , J) *COS( THETA) 

ENDIF 

ENDIF 

IF(R(I ,J) -EO.O.O)  R(I,J)*l.E-8 
20  CONTINUE 
10  CONTINUE 
RETURN 
END 
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SUBROUTINE  METRIC <  X , Y , RE , IMET , DX , DY , DYXI , DYB , DYC , DYF , 

&  ETAX , DYRE , IMAX , JMAX, JHM1 ) 

C  CALCULATES  DERIVATIVES  XI  WRT  X,  ETA  WRT  Y,  AND  ETA  WRT  X 
C  XI  WRT  X  =  l/DX.ETA  WRT  Y=1/DY,AND  ETA  WRT  X=  -DYXI/ <DX«DY) 
DIMENSION  X ( IMAX) ,Y<IMAX, JMAX) , DX(IMAX) ,DY( IMAX, JMAX) 
DIMENSION  DYXKIMAX, JMAX)  ,DYB(IMAX, JMAX) 

DIMENSION  DYF < IMAX , JMAX ) , ETAX ( IMAX , JMAX ) , DYRE < IMAX , JMAX) 

C 

DX(1)“.5»(-3.«X<1)*4.«X(2>-X(3)) 

DO  10  1*2, IMAX 
IM1*I-1 
IM2*I-2 
IP*I*1 

IF (I. EO. IMAX) THEN 

DX(I)*.5»(3.»X(I)-4.»X<IM1)+X(IM2)) 

ELSE 

DX(I)*(X<IP)-X(lMl))/2. 

ENDIF 

DYXI(I,1)=0.0 
IF  (I.GT.IMET)  THEN 

C  DYXI < I , JMAX) =Y  < I , JMAX) -Y( IM1 , JMAX) 

DYXI ( I , JMAX )*.5»<3.»Y(I, JMAX) -4.*Y(IM1, JMAX ) ♦ Y ( IM2 , JMAX ) ) 
ELSE 

DYXI ( I , JMAX) * Y  < I , JMAX ) • .5»DX ( I ) /X ( I > 

ENDIF 

DO  10  J*2, JMM1 
JP=J*1 
JM1-J-1 
JM2= J-2 

DY(I,J)*(Y(I,JP)-Y<I,JMl))/2 
IF (I.GT.IMET) THEN 
C  DYXI (I , J) *Y(I, J) -Y(IM1 ,  J) 

DYXI  (I ,  J)  =  .5*  <3.  »Y(I,  J)  -4.  »Y(IM1,  J)  ♦Y(IM2,  J) ) 

ELSE 

DYXI(I,J)=Y(I,J)».5»DX<I)/X<I) 

ENDIF 

DYB( I,J)-Y(I,J)-Y(I,JM1) 

DYC=(Y(I,  J°)'-Y(I,JM1) )  /2. 

DYF (I,J)*Y<I,JP)-Y(I,J) 

ETAX(I,J)=- DYXI (!,J)/(DX(I)*DY(I,J)) 

10  DYRE (I ,  J)  =  l ./ <DYC»RE) 

RETURN 

END 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  ERROR  IN  U-UEXACT 
C 

SUBROUTINE  ERROR <U,UEXACT, IMAX, JMAX) 

DIMENSION  U < IMAX, 0 : JMAX* 1) , UEXACTC IMAX, 0: JMAX* 1) 

URMS=0.0 
DO  10  1*1, IMAX 
DO  20  J  =  1 , JMAX 
UER*U<I,J)-UEXACT(I,J) 

URMS=URMS*UER»*2. 

20  CONTINUE 
10  CONTINUE 

RMSU=SQRT (URMS/ ( IMAX* JMAX) ) 

WRITE<6, » )  'RMS  VALUE  OF  U-VELOCITY  ERROR  =',RMSU 

RETURN 

END 
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SUBROUTINE  NWP ( U , DUUU , H , IMAX , JMAX , ICOL , WHO , WH2 , WH3 ) 

DIMENSION  U<IMAX,0:JMAX1)  DUUU(JMAX) 

RMSDU*  0.0 

RHSDUU=0.0 

RMSDUUU=0.0 

H0=0.0 

H2=0.0 

H3=0.0 

WRITE  <6,100) 

100  FORMAT  <3X, 'ETA', 07X, 'DU' ,9X, 'DUU' ,10X, 'DUUU', 7X, 'DUUU**2' ) 
DO  10  1*1, IMAX 
C  I*ICOL 

IF ( I .EQ.2.0R. I .£0. IMAX)  WRITE(6,*) 'I*'  ,1 
DO  20  J=l, JMAX 

IF  (J.EQ. 1 .OR. J .EQ. JMAX)  THEN 
IF  U.EQ.l)  THEN 

DU=.5«(-3.»U(I,1)*4.«U(I,2)-U<I,3)) 

DUU=U(I , 1) -2. *U (I,2)*U(I,3) 

ELSE 

DU* .5» (3. »U (I , JMAX) -4 .  »U  (  J , JMAX-1 ) ♦  Uf  I , JMAX-2) ) 

DUU*U(I, JMAX) -2.  »U(I ,  JMAX-1)  ♦Ud ,  JMAX-2) 

END  IF 
ELSE 

DU=.5*<U<I,J+1)-U<I,J-1) ) 

DUU=U(I,J*1)-2.»U<I,J)»U<I,J-1) 

END  IF 

IF  <  J . LE . 2 . OR . J . GE . JMAX- 1 )  THEN 
IF  (J.LE.2)  THEN 

DUUU< j)=l5« (-3. »U(I,J*4)*14.*U(I,J»3) -24. «U(I , J+2) 

1  ♦18.»U(I,J+1)-5.»U(I,J)) 

ELSE 

DUUU( J) * .5* (5. »U(I , J) -18. «U(I , J-l) ♦24. »U(I , J-2) 

1  -14.«U<I,J-3)*3:»U(I,J-4)) 

ENDIF 

ELSE 

DUUU< J) *.5»<U<I,J+2)-2.»U<I,J+l)+2.«U(I,J-l)-U(I,J-2)) 

ENDIF 

H3=H3*DUUU( J) «»2 
H2*  H2*DUU«*2. 

HO*  HO  ♦  DU»*2. 

RMSDU*  RMSDU  ♦  DU»*2 
RMSDUU*  RMSDUU  ♦  DUU»«2. 

RMSDUUU=RMSDUUU+DUUU< J) »»2 
IF( I .EQ.2.0R. I .EQ. IMAX) 

1  WHITE  <6,200)  J,DU.DUU,DUUU<J) ,DUUU(J)**2. 

200  FORMAT  <1X.I5,4E13.5) 

20  CONTINUE 
10  CONTINUE 

H  =  WH3-H3  ♦  WH2-H2  ♦  WHO»HO 
RMSH=SQRT (H/ < IMAX* JMAX) ) 

RMSDU=SQRT (RMSDU/ <IMAX» JMAX)) 

RMSDUU=SQRT<RMSDUU/ < IMAX* JMAX) ) 

RMSDUUU*SQRT(RMSDUUU/ < IMAX»JMAX) ) 

WRITE (6 , 300 )  H, RMSH , RMSDU , RMSDUU , RMSDUUU 
300  FORMATUX,  'THE  H=' ,E15.9,3X, 'THE  RMS  OF  H=',E13.5,/, 

1  IX, 'THE  RMS  OF  DU*' ,E13.5,3X, 'THE  RMS  OF  DUU*',E13.5, 

2  /,1X, 'THE  RMS  OF  DUUU*' ,E13,5,/> 

RETURN 

END 
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C  THIS  SUBROUTINE  CALCULATES  THE  POHLHAUSEN  APPROXIMATE  SOLUTION 
C  FOR  INCOMPRESSIBLE.  FLAT  PLATE. 

C 

SUBROUTINE  POHL ( Y , DELTA , UE , IMAX , JMAX . UEXACT ) 

DIMENSION  Y<IMAX,JMAX),DELTA<IMAX>,UEXACT<IMAX,0:JMAX*1) 

AMBDA=0 .0 
A*0. 0 

B=2.0*AMBDA/6. 

C=-AMBDA/2. 

D=-2. +AM8DA/2. 

E=l-AMBDA/6. 

DO  30  1*1, IMAX 
UEXACT(I.O) *0.0 
30  CONTINUE 

DO  10  1=2, INAX 
DO  20  J=1 , JMAX 
IF(Y(I, J) .LE.DELTA(I) )THEN 

UEXACT ( I ,J)=<A+B»<Y<I,J) /DELTA  ( I )  ),+C»  <  ( Y  ( I , J) /DELTA( I ) ) «*2. )* 
1  D*( (Y(I, J) /DELTA( I ))»»3.)-*E»((Y<I,J) /DELTA ( I) ) *«4. ) ) 

ELSE 

UEXACT  < I , J) =UE/UE 
END  IF 

20  CONTINUE 
10  CONTINUE 

DO  60  J=1 , JMAX 
UEXACT  <1 , J ) * UEX ACT ( 2 , J  > 

60  CONTINUE 
RETURN 
END 
C 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  NEW  GRID  GENERATION  CONTROL 
C  FUNCTION  P.  IT  IS  BASED  ON  A  GLOBAL  TRUNCATION  ERROR 
C  ANALYSIS  FOR  THE  FLOW  SOLUTION  IN  THE  TRANSFORMED  PLANE 
C 

SUBROUTINE  NEWP1  <U,P.UPLUSDC,UMINDC,DUUU,CPLUSDC, 

1  CMINDC . C , RMSDUUU ,H,K, PERCENT , IMAX , JMAX , 

2  ABAR2.NMAX.X,Y.AAA.CCC.DDD,GGG,WWW.YNIN,YMAX,IC0L, 

3  ITERATE.  UE.  DELTA.  P3,Y3.D3U,U3,  A. ZI.IR,  WHO,  WH2.WH3) 

DIMENSION  DUUU(JMAX) ,P<JHAX) 

DIMENSION  UPLUSDC ( IMAX , 0: JMAX+1 ) , UMINDC ( IMAX.O: JMAX*1) 

DIMENSION  U(IMAX,0: JMAX»i) , U3( IMAX, 0: JMAX* 1) 

DIMENSION  P3< JMAX) ,Y3(1MAX. JMAX) ,D3U(JMAX) 

DIMENSION  X(IMAX) ,Y<IMAX, JMAX) .DELTA ( IMAX) 

DIMENSION  ABARP(NMAX) ,  A(NMAX)  ,ZI  (NMAX.NMAX)  ,WWW(JMAX) 

DIMENSION  AAA(JMAX) ,CCC( JMAX) ,DDD< JMAX) .GGG(JMAX) 

WRITE(6, 100) 

100  FORMAT </.15X,' . NEW  RLAMBA  CALCULATION  ***■****«*',//, 

1  3X,'ETA',5X, 'DU' ,11X, 'DUU' ,10X, 

2  'POLD' , 10X, ' DUUU' ,7X, 'DUUU»»2' ) 

RMSDU=0.0 

RMSDUU=0 .0 

RMSDUUU=0.0 

H2=0.0 

H3=0.0 

HP2=0.0 

HP3=0, 0 

HM2=0.0 

HM3=0.0 

HP0=0.0 

HMO=0.0 

H0=0.0 

COLD=C 
DO  20  1  =  1, IMAX 
C  I=ICOL 

IF  < I . EQ . 2 . OR . I . EQ . I MAX )  WRITE(6,«)'I=  ',1 

DO  10  J=l, JMAX 

IF  ( J . EQ . 1 . OR . J . EQ . JMAX )  THEN 
IF  .(J.EQ.l)  THEN 

DU=.5*(-3.»U(I,1)*4.»U(I,2)-U(I,3)) 
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DUP= . 5* < -3 . "UPLUSDC ( I . 1 ) "4 . "UPLUSDC (1,2) -UPLUSDC (1.3)) 

DUM* ,5» ( -3. "UMINDCd , 1 ) *4. "UMINDC (1 ,2)  -UMINDCd , 3) ) 

DUU=U( I ,l)-2. *U(I ,2) ♦U (1,3) 

DUUP  =  UPLUSDC ( 1 . 1 ) -2 . "UPLUSDC (1,2) "UPLUSDC (1.3) 

DUUM=UKINDC(I , 1) -2. "UMINDC (I , 2) +UMINDC (1,3) 

ELSE 

DU=.5"(3."U(I , JMAX) -4 .  »U (I , JMAX-1 ) *U (I .JMAX-2) ) 

DUP= . 5* ( 3 . "UPLUSDC ( I , JMAX) -4 . "UPLUSDC ( I , JMAX-1 ) ♦ 

1  UPLOSDCd, JMAX-2) ) 

DUM=.S*  <3.  "UMINDCd,  JMAX) -4.  "UMINDCd,  JMAX-1)  ♦ 

1  UMINDCd, JMAX-2) ) 

DUU=U(I,  JMAX) -2.  *U(  I,  JMAX-1)  "Ud,  JMAX-2) 

DUUP= UPLUSDC ( I , JMAX ) - 2 . "UPLUSDC (I , JM AX- 1 ) ♦ 

1  UPLUSDCd, JMAX-2) 

DUUM=UMINDC(I . JMAX) -2 . "UMINDC( I .JMAX-1 ) ♦ 

1  UMINDCd.  JMAX-2) 

END  IF 
ELSE 

DU*.5*(U(I,J+l)-Ud,J-l)) 

DUP= . 5" ( UPLUSDC (I, J*l) -UPLUSDC (I . J- 1 ) ) 

DUM=  .5"  (UMINDCd ,  J*l)  -!JMINDC(  I ,  J-l ) ) 

DUU=U(I, J*l)-2. •U(I,J)*U(I,J-1) 

DUUP=UPLUSDC (I,J*l)-2. "UPLUSDC (I . J ) ♦UPLUSDC (I , J- 1 ) 
DUUM=UMINDC(I,J»1)  -  2.  "UMINDCd,  J)  ♦UMINDCd  ,J-1) 

END  IF 

IF  ( J.LE.2.0R. J.GE. JMAX-1)  THEN 
IF  ( J.LE.2)  THEH 

DUUU( J) * ,5*(-3. *U(I , J*4)*14. "U(I, J»3)-24. "U(I , J"2) ♦ 

1  18.»U(I,J*1)-5.»U(I,J> ) 

DUUUP=  .5"  ( -3.  "UPLUSDCd  ,J«4)  "14 .  "UPLUSDCd  ,  J+3)  - 
1  24 . "UPLUSDC (I, J+2) ♦18. "UPLUSDC ( I , J*1 > -5 . "UPLUSDC (I , J) ) 

DUUUM=.  5"  (-3.  "UMINDCd,  J"4)  "14.  "UMINDCd,  J"3>- 
1  24. "UMINDCd,  J»2)d8. "UMINDCd,  J-D-5. "UMINDCd, J)) 

ELSE 

DUUU(J)*.5*(5."U(I,J)-18.*Ud,J-l)*24."U(I.J-2)-14."U(I,J-3)» 

1  3. "U(I , J-4) ) 

DUUUP=. 5- (5.  "UPLUSDCd,  J)-18. "UPLUSDCd, J-l)  *24. "UPLUSDCd .  J-2)- 
1  14. "UPLUSDCd,  J-3)*3. "UPLUSDCd, J-4)) 

DUUUM=. 5"  (5.  "UMINDCd,  J) -18.  "UMINDCd,  J-l)  "24.  "UMINDCd,  J-2)  - 
1  14. "UMINDCd.  J-3)"3. "UMINDCd,  J-4)  > 

END  IF 
ELSE 

DUUU(  J)  =  .5"  (U(I,J»2)-2.*U(I,J"l)»2."Ud,J-l)-U(I,J-2)) 

DUUUP=.  5"  (UPLUSDC(  I.  J»2) -2.  "UPLUSDCd,  J*l)  *2.  "UPLUSDCd,  J-l)  - 
1  UPLUSDCd.  J-2)) 

DUUUM=.  5"  (UMINDCd,  J"2) -2.  "UMINDCd,  J*l)  *2.  "UMINDCd,  J-l) - 
1  UMINDCd.  J-2)) 

END  IF 

IF(I .EQ.2.0R. I .EQ . IMAX)  THEN 

WRITE  (6,200)  J,DU,DUU,P(J),DUUU(J),DUUU(J)""2. 

200  FORMAT  (I5.5E13.5) 

ENDIF 

RMSDUUU=RMSDUUU"DUUU(J)"*2 
RMSDU = RMSDU ♦ DU  »  »  2 
RMSDUU=RMSDUU*DUU"*2 . 

H0=H0»DU""2. 

H2=H2*DUU»»2. 

H3=H3»DUUU( J) »»2 
HP0-HP0»DUP"»2. 

HP2=HP2»DUUP*"2. 

HP3=HP3"DUUUP»»2 . 

HM0=HM0*DUM"»2. 

HM2=HM2*DUUM»"2. 

HM3=HM3*DUUUM»»2 
10  CONTINUE 
20  CONTINUE 

H=WH3»H3  ♦  WH2*H2  ♦WHO'HO 
HPLUSDC=WH3»HP3  ♦WH2*HP2  ♦WHO»HPO 
HMINDC=WH3"HM3»WH2"HM2  ♦  WHO"HMO 
DET= (CPLUSDC»»2) • (C-CMINDC) -CPLUSDC* (C*C-CMINDC»"2) * 
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X  C»C*CMINDC- (CMINDC*  »2)  *C 

RNUM= (CPLUSDC»*2) * (H-HMINDC) -HPLUSDC* (C*C-CMINDC»<*2) ♦ 

1  (C»«2)«HMINDC- (CMINDC»»2) »H 

DEN=HPLUSDC* (C-CMINDO) -CPLUSDC* (H-HMINDC) ♦H«CMINDC- 
1  HMINDC*C 

IF  (DET.NE.O.O)  THEN 
C=-.5*RNUM/DEN 
IF  (DEN/DET.LT.O.O)  THEN 
HMIN=MIN (HMINDC, H , HPLUSDC) 

IF  (HMINDC.EQ.HMIN)  C=CMINDC 
IF  ( HPLUSDC. EQ.HM IN )C=CPLUSDC 
IF(n.EQ.HMIN)  C=COLD 
ELSE  . 

DO  40  JJ=1,NMAX 

ABAR2( JJ)  =A( JJ) ♦  C  *  2KJJ.IR) 

40  '  CONTINUE 

CALL  ?CAL(ABAR2,NMAX,P3, JMAX) 

CALL  EXPGRD < X , Y3 , P3 . AAA . CCC , DDD , GGG , WWW , YM IN , 

1  YMAX , IMAX, JMAX , ICOL) 

CALL  BLIMP ( X , Y3 , IMAX , JM AX . U3 .ITERATE , UE , DELTA , RE) 

CALL  NWPCJ3.D3U.H4, IMAX, JMAX, ICOL. WHO. WH2.WH3) 
HMIN=NIN(HNINDC, H, HPLUSDC. H4> 

IF < H4.EQ. HMIN. AND. H4.NE.H)  GO  TO  220 
IF  (HMINDC.EC.HMIN)  C=CMINDC 
IF  (HPLUSDC. EQ.HMIN)C=CPLUSDC 
IF(H.EQ.HMIN)  C=COLD 
220  CONTINUE 
END  IF 
ELSE 

HMIN=MIN(HNINDC,H, HPLUSDC) 

IF  (HMINDC.EQ.HMIN)  C=CMINDC 
IF  (HPLUSDC. EQ. HMIN)  C=CPLUSDC 
IF  (H.EQ.HMIN)  C=COLD 
END  IF 

EPSILON* .005 

IF  (ABS(C) .LT. .01)  EPSILON* .01 
DC*PERCENT*C-*EPSILON 
CPLUSDC=C»DC 
CMINDC=C-DC 

WRITE  (6,300)  H, HMIN, HPLUSDC, HMINDC.C.K 
300  FORMAT  ( IX . ' H= ' , E15. 9, / , 

1  lX,'HMIN=',E15.9,y, 

2  IX, ' HPL=' ,E15 . 9, / , 

3  lX.'HMN*' ,E15.9,/, 

4  IX, 'NEW  RLAMBA*' , E15.9, / , 

5  IX. 'ITERATION  NUMBER,  K=',I5,/> 

RMSDUUU* SORT (RMSDUUU/ (IMAX*JMAX) ) 

RMSDUU*  SQRT(RMSDUU/(IMAX»JMAX) ) 

RMSDU*  SQRT(RMSDU/ (IMAX*JMAX) ) 

WRITE  (6,400)  RMSDU . RMSDUU , RMSDUUU , COLD 
400  FORMAT  (IX, 'RMSDU  =  ' ,E13. 5, 3X , 'RMSDUU  *' .E13.5.3X,  'RMSDUUU 
1  ,E13.S,/,1X, 'RLAMBA  OLD* ' , E15 . 9, / ) 

RETURN 

END 
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l J  Hr  vs  v 

This  e-t-u d y,_d e  v  e  1  o  p  s  sn  adapt  ive— gr  id  method  which  mini¬ 
mixes  the  trnncation  error  in  the  finite-difference  solution. 
The  study  solves  compressible,  steady-estate,  boundary-layer 
equations  assuming  pe r f e c t—g a s  flow  over  an  isothermal  wall. 
The  Do r odp i t syn ,  compressibility  transformation  change?  the 
boundary-layer  equations,  as  expressed  in  two-dimensional, 
cartesian  coordinates,  into  an  incompressible  form.  The 
equations  are  theh  transformed  into  variables  of  a  compu¬ 
tational  plane.  Implicit  Sue  c  e  s  s  ive—Ove  t—'fce  1  axa  t  ion  (SOR) 
solves  the  f in i t e— d i f f e r enc e d ,  computational,  boundary— layer 
equations.  Comparison  of  the  computed  solution  for  incom¬ 
pressible  flow  over  a  flat  plate  to  Blasius',  exact  solution 


shows  the  boundary-flayer  code 


accurate. 


The  adaptive-grid  method  uses  Powell's  method  to  opti¬ 
mize  the  solution  grid  by  minimizing  the  sum  of  the  third 
derivative  in  the  computational  plane  of  the  tangential  velo¬ 
city  component.  Powell's  method  finds  the  grid,  control 
function,  Q,  in  anXelliptic,  grid  equation,  Yt  (  +  QY^-fO, 
which  minimizes  a  specified  function.  The  grii  equation 
generates  the  grid  spacing  at  the  end  of  the  plate.  This 
spacing  is  then  strehmwise  scaled  across  the  remaining  grid. 
Minimizing  the  sum  of\the  square  of  the  third  derivative  in 
the  computational  plane\of  the  tangential  velocity  component, 
U^.  ,  over  the  entire  d<im<in  decreases  the  truncation  error 

the  best  of  the  functions  t  e  s  Ve  d  t-v  This  study  tests  the  sum 
of  the  squares  of  the  first  second,  and  third  derivative  of 
the  tangential  velocity  as  minimized  functions.  The  accuracy 
of  the  optimized,  adaptiv e— g rid  solution  is  greater  than  the 
original,  fixed— grid  solution.  The  study  then  applies  the 
optimization  to  supersonic  and  hypersonic  problems.  The 
computed,  ad ap t i ve— g r id  solutions  show  good  correlation  with  , 
theoretical  models  for  supersonic  and  hypersonic  flow,  devel¬ 
oped  by  Van  Driest  and  Eckert. 
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